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I. INTRODUCTION AND BACKGROUND 


A. BACKGROUND 


Suppose that an underwater acoustical source emits a 
Signal at a known time. Suppose that the times of arrival 
of that signal at the hydrophones ona three dimensional 
sensing array are also known. Finally suppose that the 
speed of scund in water is modelled as being homogeneous 
over time and horizontal displacement, varying only with 
depth. Ther the angle of elevation (å), and the time (T) of 
the arrival of the signal at the acoustic center of the 
hydrophone array can be determined. If the :elationship 
between the speed of sound and tke depth under water is 
known exactly, then the angle A and the time T can be used 
to trace the signal trajectory back over its ray pata 
(Called ray tracing) to determine the original position oZ 
the acoustical source [Ref. 1]. However, full realization 
of this method in actual hydrophonic tracking is prevented 
by two primary sources of inaccuracies in the process. The 
first and probably greatest problem is that tie speed of 
sound profile can be approximated at only a few locations, 
usually at great cost to achieve even moderate accuracy. In 
addition the profile is certain to flucuate over time and 
location: The second problen, confounding the first, is 
that there may be innacuracies, of uaknown size, in the time 


data values recorded by the sensing array. 


B. PURPOSE 


As noted in [Ref. 1] the procedure of determining a 
Sound source position by ray tracing is very sensitive to 


even small errors in the angle of elevation or speed of 


10 


sound at the sensing array. Of course ray tra- ing is also 
highly dependent on the accurate determination of the the 
transit time from source to array. The purpose of this 
study is to develop an appropriate model of the timing data 
observed in the hydrophonic tracking problem. Tie objective 
or the desired model is to eStimate the error in the time 
values, and produce improved estimates of the initial angle 
and ray transit time, so as to reduce the effect of those 
errors on the ray tracing procedure. 

In pursuit of these objectives this thesis first reviews 
the currently used models, which are called the NAVY and 
NAVY_A methods. Then four alternative models are developed, 
cece Mice, aoe Ce, Nel. P. and Moias. methois. Finally 
the performance of all six models are evaluated and compared 
usiny simulation studies. The comparisons are made under 
the idealized conditions of a known linear spesd of sound 


versus depth relationship. 


C. TRACKING RANGE CONFIGURATION 


The tracking range which supplied data for this study 
consists of several separate three dimensional hydrophone 
arrays sitting on the sea botton. They are lared out ina 
rough grid, each array being approximately 750) feet fron 
the next, so that the sound source being tracked is never 
more than about 5000 feet from the nearest array. The 
arrays are at depths cf roughly 1000 to 1300 feet. 

Each array (see figure 1.1) has four indeperdent hydro- 
phones defining an orthogonal coordinate system. The phones 
are referred to as the x,y,z, and c hydrophones, and are on 
Tour adjacent vertices of a cube (see figure 1.2 = with sides 
of length D (usually 30 feet). The arrays are linked to 
shore based computers by electronic cable. IO nO E 


the lccal coordinate system oí each array is at the acoustic 
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Figure 1.1 Acoustic Signal Detection by 
Three Dimensional Hydrophonic Array. 
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Figure 1.2 Geometry of Hydrophonic Arriys. 


The sound source is equipped with a clock synchronized 
with the shcre based computers and emits a signal at speci- 
fied intervals. The times of receipt of those signals at 
the four hydrophones are recorded, and the corresponding 
travel times are calculated by subtraction of the signal 


emission time. 


D. APPARENT POSITIONS 


The first step in estimating the position of a sound 
Source is to do so under the assumption that th» sound wave 
travelled its entire trajectory through water which had a 
constant speed of sound. The result, called the apparent 


position, is obviously erroneous, but is then corrected by 


D 


the ray tracing procedure (a description of waren omer 
later)... The constant velocity value used is usually the 
estimated speed of sound V at the depth of "he sensing 


array. 
The constant velocity assumption implies that the wave- 


front is an expanding sphere. Therefore the apdarent posi- 
tion is calculated using Simple spherical ecuations 
involving squared distance calculations. Specifically, is 


Tx is the travel time recorded for the x-phone, then  VeTx 
is the distance between the apparent position and the 
x-phone. This distance is also egual to the usual geometric 
distance between the two positions 
( X EY E) (apparent position) 

and (D ,-D ,-D) /2 (x hydrophone position). 
Equating these two squared distances, and equating their 
counterparts for the other three hydrophones, tie equations 


(1.1) are oktained. 


(x= a2) + (v4 7/2)? + 0240/2)? = ve N; 
(A A A ene I 
(X+ Dae ey + D/2 L +o ZoD je = P 
( X * D/2 E + ( Y + D/2 E + ( Z + D/2 y = cane 


Assuming that the times Tx, Ty, Tz, and IT: are known, 
(1.1) is a system of four equations in three unknowns X, Y, 
and Z. This overdeternined system will, in genecal, have no 
exact solution. In fact, even 1f the time values were 
exactly correct, the system would still have no exact solu- 
tion: This is because the equations corresoná to the 
straight line ray paths due to the constant velocity assunp- 
tion, whereas the time values correspond to tie true ray 
paths which are not straight due to the “aetual vara c s 
velocity of sound alcng the ray path. This is a subtle, but 


very imrortant point. 
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DegirorowfoutMone Of the equations arbitrarily, so as to 
reduce the system to three equations in three unknowns, is 
to throw away information from one of the hydrophones. The 
psuedo solution currently utilized is to subtract the fourth 
equation from each of the first three, yielding a system of 
three equations in three unknowns which allows an exact 
soluticn involving information from all four hydrophones. 
However that solution will not, in general, satisfy any of 
the original four equations, and is only one of many reaon- 


able ways to choose an approximate solution. 


E. INITIAL ELEVATION ANGLE AND RAY TRANSIT TIME 


Assuming that a solution (Xa,Ya,Za) has bee. determined 
for the apparent position, then the initial angle of eleva- 


tion is just the angle of elevation of that solution, given 


ba (1.2). 
A, = arcsin Z J> | x + y? + z^ Ç le) 
1 a a a a 


The objective is to find an apparent positio1 (Xa,Ya,Za) 
such that (1.2) computes an angle which  approximates the 
physically correct elevation angle as closely as possible. 
Me solution (Xa,Ya, Za) and the times Ix, Ty, Iz and Tc can 
then ke used to determine an appropriate value for T, which 
is the ‘ray transit time', or time of arrival of the sound 
Wave at the acoustic center of the sensing acray. The 
Currently employed method uses the proportional relationship 
or ecudation (1.3), where R and Rc are the ranjes from the 
apparent position to the acoustic center and to the c hydro- 


phone respectively, as in (1.4). 
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R T_ R (AMEN 








e R C 
= V = —  — — = 
c C 
Z 
R = X t Y + Z 
a 


(1.9) 


2 2 2 
(X *D/2) è (Y +D/2) + (Z_+D/2) 


ye) 
II 


F. RAY TRACING 


Whichever method is selected to produce tie apparent 
position (Xa,Ya,Za), it is transformed into the estimate of 
the true sound source position (X,Y,Z) by the procedure of 
ray tracing? When there is velocity layering ir the water, 
the ray path is no lcnger a straight line. This is treated 
using repeated applications of Snell's Law (Rei. 2 p.131], 
starting with the layer of water in which the array sits, 
and backtracking upwards through successive velo: ity layers, 
until the estimated Fay transit merr Prs Con Omen 

The layering effect is artificially induced by the limi- 
tation that the speed of sound can be estimated at only a 
finite number of depths, the result of which is commonly 
called the water column. For example, at the tracking range 
studied the speed of sound is measured every 25 feet, 
starting at the depth of 12.5 feet. Hence, for example, the 
third layer from 50 to 75 feet deep is assumel to havea 
constant speed of sound equal to that measured at 62.5 feet. 

The first layer processed [Ref. 3 p.4] is the partial 
layer lying between the array and the deepest layer boundary 
that is shallower than the array, with thickness 2Z1 (see 
Liguria 1.35). 

The incremental slant range in the first Layer is S1 
given by (1.5), where Ai is the initial elevation angle 


estimate. 
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Figure 1.3 Ray Tracing an Apparent Position. 


SH Sin (a,) 


travel time 


(1.5) 


The incremental in the first Layer is T1 


given by (1.6), 


where Y1 is 


the velocity estimated for the 


layer in which the array is situated. 


aer 


T) = S)/v 


1 (1.6) 


The incremental horizontal distance travellel by the ray 
in the first layer is H1! jiven by (1.7). 
Hi = S; cos (A) ebe 


To determine the angle of elevation in the next layer, 
Snell's Law (1.8) is applied, where V2 is the speed cf sound 


estimated for that layer. 


cos (A. ) eos (a) 
1 E ee ia 
Vi Vo 
When (1.8) is sclved for the cosine of tie angle of 


entry into the next layer, (1.9) is obtained. 


| e V. cos(A,) 
COS A») = ere (1 9) 


E 


The procedure of computing the incremental values of 
Slant range, time and horizontal distance are now repeated 
for the second layer. The overall procedure is repeated 
upwards through layers 2,...,n , where n is the first layer 
in which the sum of the incremental travel times exceeds the 


total ray transit time, as in (1.10). 


(ës Sab e EE T GL TS 


in the last, uppermost layer the values Tn, Sn, Hn, and 
Zn must ke adjusted to compensate for overshooting the total 
time T. The values Hi and Zi (130%. are then accumu-= 
lated to get (1.11). 


n 
H = La H; z = E e (ap 
1= n= 


Now the raytraced position estimate is given by (1.12). 
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ZAS a [SETA 
S a a (1. 12) 

Y = co Vd y? 

a 


= 

The sensing array is usually not aligned with the coor- 
dinate system of the overall tracking range. Therefore the 
apparent position, which is in terms of the local array 
coordinates, must be changed by a suitable geomstric trans- 
formation prior to ray tracing SOMA Sto account for the 
angle of tilt at which the array sits on the sea bottom. 
meter ray tracing, the position estimate must be again 
transformed to account for rotation of the array about its Z 
axis away from a position which is aligned with the range 
coordinate axes. Finally a simple translation must be 
applied to reference the position estimate to the range 
coordinate system origin vice the acoustic ceter of the 
array. The end result is a position in terms of the overall 
range coordinate system. These transformations are not 
given here because they are used after the estimation of the 
initial angle and time, and hence do not affect the accuracy 
of those estimates. See [Ref. 3] for further details. 


G. DISCUSSION 


If the velocity versus depth relationship is smooth and 
estimated accurately, then the ray tracing procedure is 
Surprisingly robust with respect to the thickiess of the 
layers. For example, if velocity is linear versus depth, 
and is known exactly, then the exact hydrophone times, ray 
transit time and initial angles can be computed [Ref. 4] for 
any given sound source position. Then tke ray tracing 
procedure, with layers as thick as 25 feet and targets as 
far away as 3000 feet, Still estimates positions to within: 
inches of each of the true coordinate values. This seems to 


indicate that errors resulting from position estimation are 
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not due to the approximation by layers, except possibly when 
there are radical changes in the velocity pattern within 
single layers such as frequently occur in laye:s near the 
water surface. Rather such errors apparently ars: due mostly 
to inaccuracies either in the estimation of the speed oz 
sound profile itself, or in the initial anjle ani rransms 
tine estimates. This study will focus on those 2rrors which 
are involved in the time values observed at the nydrophones, 
and attempt to produce methods of initial angle and transit 
tine estimation which reduce the effects of those errors on 


the overall position estimation procedure. 


20 


wee BASIC CONCEPTS 


The method currently used for estimation of an initial 
angie and ray path transit time focuses on the four spher- 
ical equations (1.1) given in Chapter I. As previously 
discussed, these equations have no exact solution because: 

1. they form an overdetermined system of four equations 
in three unknowns; 

2. the time values recorded at the hydrophones € may be 
innacurate due to unknown sources of error in the 
observation process; and 

3. even if the time values were exact, they correspond 
to a nonconstant velocity profile, and so will not be 
be correct for the constant velocity geometry (spher- 
ical wavefront) assumed by the equations. 

As previously mentioned, the psuedo solutidn chosen by 
the current method is to subtract the fourth spherical equa- 
tion from each of the first three, and solve the resulting 
system of three equations in three unknowrs. This rethod 
has the beneficial quality that information is r2tained fron 
all four hyarophones, whereas to just drop the fourth eyua- 
tion (or any one of the equations) without the initial 
Subtraction would cause complete loss of the irformation 
from the data recorded by one of the hydzophones. However 
it is important to realize that the solution thus obtained 
does not actually satisfy any of the originai Tour 
equations. 

It should be noted that the development of tais solution 
in [Ref. 3] is done entirely from a geometrical point of 


view, and does not mention the system of foir spherical 
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eguations. The text of [Ref. 3] does not drawn Cren EOTS 
the fact that tne sclution developed 1S just one of many 
plausible choices, none of which will satisfy all “our 
spherical ccnstraints. Therefore the solutior chosen is 
treated as though it were the exact solution, only subject 
to errors in the observed hydrophone time values. However, 
even with exactly correct time values, this currently 


enployed method will not yield the true elevati»n angle and 


ray transit tine. This is due to the assumption of a 
constant velocity vice the true nonconstant velocity 
profile. This conflict introduces an automatic bias in the 


initial angle and time estimates currently used for ray 


tracing: 


B. COMPUTATIONS 


To simpiify notation, let (X,Y,2) be the codrdinates of 
the apparent position which was formerly denoted (Xa,Ya,Za). 
Then when the current solution is applied, the first step is 
to subtract the fourth equation from eacn of the other 


three, which produces the equations (2.1). 


(Cx DA? ys - ( X + D/2 ys = v* ( r? - T^ ) 

5 I a x LZ age 
(Y - D/2 ) - (vY+ D/2 ). = v* ( SS - T. ) 
(2-0/2 )°- = tZ 4 B2) PDP R 1s 2 it 


The solution to these are easily obtained, as in (2.2). 
2 2 E 


X = N O NE, 
a. Se (232) 
Y = y ( T^ E ) 2 
e y di D 
2 2 2 
Z = — fm 
Wu n os 


Then the initial elevation angle estimate is (2.3). 


Fa | me 2 2 
=" Warcsoan (=/ EC ) (3) 
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Mecca pa h transit time to the acousti> center is 


(2.4), where R and Rc are as defined in Chapter I by (71.4). 


Die E R (2.4) 


This method of determining the apparent position shall 
hereafter be referred to as the "Navy Method', or 'NAVY' for 


short. 


C. ADJUSTMENTS TO THE ORIGINAL SOLUTION 


Experience has shown that the NAVY method produces an 
apparent position estimate which usualiy can be improved by 
an adjustment which is described in this section., 

The cosine of the angle between the i-th axis and the 
straight line fron the origin out to the apparent position 
ìs called the i-th direction cosine Ci. Itis a fact of 
geometry that the sum of the sguares of the three direction 
cosines must equal unity. Therefore the method is now 
adjusted to reflect that constraint. 

The direction cosines used are the,angles nade by the 
ray path at the c hydrophone. Therefore the (X,Y,Z) coordi- 
nates calculated by the original method are first translated 
to coordinates referenced to the c-phone as the temporary 


SEÑA na, aS in (2.5): 
X. = X + D/2 Y, - Y + D72 EA + D/2 (2.5) 


Therefore the three direction cosines arz given by 
(22-63 
E = x = = 
SARA E A VE C, Z/ VT, 


x 


(ZS) 


The denominators in (2.6) are all VeTc because that 1s 
the range £rom the apparent position to the c hydrophone, as 
estimated by the time from the c hydrophone. Ideally these 
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cosines should add to unity when squared. Tnerə Fore i£ ce 
is the ‘direction cosines correction' factor defined by 


(2.7), then DCC should be close to one. 


2 2 2 
DCC = 
c. uu Ge (2.7) 


Deviation of DCC from unity is interpreted as an indica- 
tion of receiver timing errors, array malformation or 
invalid data at one or more of the hydrophoies [Ref. 3 
poco Currently Sai Dee lies outside tae interval 
(0.98,1.02), the data is discarded as being excessively full 
of error. The direction cosines of the remaining acceptable 
data points are rescaled using (2.8) to assure satisraction 


of the direction cosines constraint. 
C = O, DEG cr = ' = m 
2 x 7 Š Si DEE c, C, / vee (2.8) 
A corrected set of new coordinates are zomputed by 
(2.9), still being referenced to the c nydrophone. 
Do CERO: E C y E cu z (259) 


These are then translated by (2.10) to coordinates 


referenced to the acoustic center. 


X= X_ - D/2 Y = Yee ee Z = Z_ - D/2 (2.10) 


This adjusted method of determining the aparent posi- 
tion shall hereafter be referred to as the "Navy Adjusted 
Method', or 'NAVY A' for short. 


D. DISCUSSION 


When a sound source is within the detection range of 
nore than one sensing array, each array produces timing 


datas Ihe data from each array can be processed to produce 
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a position estimate. Ileally these multiple estimates of 
position wiil be in reasonably close agreement. However 
experience with actual tracking data has shown that this is 
not the case in many of the multiple detection opportuni- 
ties. This tendency toward disagreement between multirle 
estimates of the same position is commonly called the cross- 
over, or crosstalk, problem. This problem often occurs when 
the sound source is moving away from the tracking domain of 
one array into the tracking domain of another. This study 
focuses on improvement of the initial angle and time esti- 
mates, which hopefully will help alleviate the crossover 
Eroblen. 

The current chcice of a 'best' compromise solution 
appears to be based on reasons of simplifying geometry and 
calculations. These are worthwhile objectives, put do not 
in themselves reflect the need to estimate accurately the 
initial elevation angle and ray transit tine. Since there 
exist physically correct values for both the aigle and the 
time, those values will produce the exact position after ray 
tracing, provided that the velocity -profile is known 
exactly. The desire then is to estimate these true values 
as accurately as possible. 

The guestion at this point is whether or not the direc- 
tion cosines adjustment causes the estimated apparent fosi- 
tion to be closer to the true apparent position. Experience 
has indicated that it does [Ref. 3 p.C-7]. However the 


effect of the adjustment can be interpreted in terms of the 


original four spherical eguations (1.1). The rescaling of 
the direction cosines given by (2.6), so as to assure that 
their squares add to unity, is eyuivalent to rescaling the 


quantities in (2.5) so as to assure that their squares add 
to (VeTc)2. That is exactly the constraint stated by the 
fourth srherical equation of (1.1). So the effect of the 
adjustment is to require that the EO UI E equation, 
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concerning the data at hydrophone c, is exacti; satisfied. 
This reguirement will, in general, assure that the other 
three equations are not satisfied. Since experience shows 
that the adjustment often improves the soiution, this seens 
to imply that the fourth equation is somehow mo: e important 
than the other three. Or it may just be that exact satis- 
facticn of ore of the equations usually assures a reasonabiy 
gool ccmpromise solution. 

Tc summarize, the NAVY method provides a useable 
apparent position suitable as input for ray tricing. But 
the direction cosines adjustment used in the NAVY_A method, 
for reasons not understood at this time, appears to improve 
tnat position as indicated by test results. Ta ose results 
are supported by the results of this thesis (see Chapter V). 
However, aS will be demonstrated by the example considered 
in the next section, the DIC correction factor ol the NAVY_A 
method has an effect which must be something more than just 


the smoothing of timing errors. 


E. A COMPUTATIONAL EXAMPLE 


For the purposes of illustration and comparison, suppose 
that a 30 foot sensing array is at a depth of 1300 feet, 
that the coordinate system origin is at the array acoustic 
center, and that the array arms are parallel to the coordi- 
nate axes. This implies that the four hydrophones are in 
the positions: 

: (lp 150 219) 
re 
: (7-157 7-15 $-— 15) 
: (-15 , -15 , -15) š 
If there is a sound source known to be in position 
( 1000 , 3900 900% 
then the depth of that source is 1300 - 900 = 40) feet. 


QO NK > 
ee 


2.6 
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Figure 2.1 


Figure 2.1 shows the estimated sound yelo: ity prorile 
which was estimated for the data used in the course of this 
thesis. As can be seen, the profile is primarily linear at 
depths greater than 100 feet. The profile at deptas greater 
than 100 feet is reasonably approximated by the linear 
relationship 

V = 4840.7 + 0.03314 e DEPTH . 

Therefore suppose that in this example probiem the 
velocity profile is known exactly, and is given dy the above 
linear relationship. 

Under these circumstances, with known linzar velocity 
profile, and known sound source location, the exact times of 
arrival of the sound wave at the four hydroptones can be 
computed using the methods set forth in Appendix A. Those 
exact times (in seconds) are: 

TX $205 026773656 595 Ty 4$ 7056731225095 
IZ 3990:675229]399» Tc : 0.079952 R 

The corresponding exact values for the initial elevation 
angle, ray transit time and resuiting apparent jcS5osition are 
also directly computable. Those computations will hereafter 
be be known as the EXACT method, and will »roduce the 
correct true position after ray tracing. The results of the 
EXACT method are given in Table I, along with the corre- 
sponding apparent position estimates produced waen the two 
methods, NAVY and NAVY_A, are applied to the (errorless) 
time values. 

At first glance the differences in Table I might seen 
rather small. However it is important to recall that these 
are produced under the ideal conditions of a very smooth and 
exactly known velocity profile. These idealizations are far 
from the realities of a nonlinear velocity profile which is 
estimated by a procedure involving errors which are unknown 
and protably significant: Such realities nigh: well cause 


the differences in Table I to be Significantly larger. The 
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| TABLE I 


Single Example Comparison of 
NAVY and NAVY_A fethods 


| 

x 

| Method n° tassa, “Š Pu dpe 

| NAVY 0.67526969 15. 26459 

| NAVY_A 0. 67527027 15. 26461 

EXACT 0. 67527043 15. 27002 

| Apparent Position Estimate 

| NAVY (edowse957., 2O0l7.874 , 868.143 ) 

| NAVY_A ( 1006.087 , 3018.259 , 868.255 ) 
EXACT ( 1005.966 , 3017.899 , 868.555 ) 


do o ec taa ROO U... U... is, A js O pts A ps io MM E, A, A rn, ADA p) 


1 am << AE I Q  — m — O “mama oy 





nature and size of those differences remain difficult to 
determine urtil more is known about the the velocity profile 
estimation errors and their effect on the position esti- 
mating process. In any event even the small differences in 
Table I might be magnified during the ray tracing process 
Mer tie conditions cf a realistic velocity profile. 

The differences illustrate the very important point that 
the direction cosines adjustment causes Changes in the esti- 
mates even when the time data is free of all erzor. Hence 
M EUuewvarMcn rrom "1.0 of the correctior factor DCC is not 
gust due to array  malfunction, receiver timlig error or 
other sources of non-valid data as previously assumed. 

In the example above, the NAVY A method produced a 
slightly better time and angle value than the JAVY method. 
However, in this same example the NAVY_A method produced an 
apparent position estimate which is slightly farther away 
from the EXACT answer than the position estimated by the 
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NAVY method. This is only one example, and its results 
should not be generalized. However it illustrat>s the point 
that apparently the true effect oí the adjustment may not be 
well understood. 

Furthermore, the adjustment seems to olace heavy 
emphasis on the time value recorded at hyirophone c. 
Therefore the effect of the adjustment may well depend 
largely on the accuracy of that one particular data value, 
which is a relatively unbalanced dependence in the presence 


of data errcrs. 
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III. PLANAR WAVEFRONT MODELS 


A. THE PLANAR WAVEFRONT ASSUMPTION 


When a sound wave travels in constant veloci: y water, it 
expands in the shape of a sphere. If the velocity profile 
is variatle instead, but is reasonably well behaved versus 
depth, then the expanding wave is a smooth distortion of a 
Spherical surface. In either case, if the wave has 
travelled a long distance when it arrives at a hydrophonic 
array, then that small piece of the wavefront #hich passes 
through the 30 foot cube spanned by the array may be approx- 
imated reasonably by a flat planar surface. This approxima- 
tion is the basis of the planar wavefront models developed 


in this chapter. 


B. PLANE EQUATIONS 


A plane in space is fully defined by ideitifying any 
Pornt {x0,70,20) on the plane, and also a vector (C1,C2,C3) 
of unit lenyth which is perpendicular to that plane. The 
vector is called the unit normal vector for that plane. 
Then any point (X,Y,Z) on that plane must satisfy the equa- 
tion of the plane, namely (3.1). 


IS e laaa, E) EO 013.1) 


The perpendicular distance from the plane t» any point 
(X1,Y1,21) not on the plane is the absolute value of (3.2). 


cnm = X. O) + C. (V 


1 1 0 > - Y. ) + G, ( 2) - Zg ) (3.2) 


I 3 1 O 
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C. BASIC MODEL EQUATIONS 


Consider a coordinate system whose origin is at hydro- 
phone c, and whose axes are aligned with the three array 
arms. Let C1, C2 and C3 be the direction cosines on the X, 
Y and Z axes respectively for the vector from the origin to 
the apparent sound scurce position. Then  (€1l,C2 ACEITES 
itself a vector, of unit length, which is perpendicular to 
the planar wavefront emanating from the soind source. 
Therefore (C1,C2,C3) can be used as the normal vector for 
the wavefront plane. 

In the coordinate system referenced to the c-rphone as 
the origin, the acoustic center has coordinates (D,D,D)/2, 
where D is the length of an array arm. When the soundwave 
plane arrives at the acoustic center, it will have the 
eauation (3.3). 

Cc, xX + C, Y + C, Z = D ( C, + C, + C, ) / 2 (3.3) 

The x-phone has coordinates (D,0,0), and he distance 
between it and the soundwave plane at the acoustic center is 
(3.4), which then simplifies to (3.5). 


C, ( D2 - D) +C U ( po PU op s C n oO KB 


2 3 ER 


D (C-C. +C +C)? 


1 2 3 (325) 


This distance is measured ina direction p2rpendicular 
to the wavefront plane, and so is measured in the direction 
of travel of the soundwave. Therefore the same distance is 
also egual to (3.6). 


V oT. -T) (3.6) 
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uo) 1S thesveloeaty of sound at the a>ray, Tx is 
the time of arrival of the sound wave at the x-phone, and T 
is the time of arrival of the sound wave at the acoustic 
center. The term "distance" is used loosely in (3.5) and 
(3.6), because these quantities nay be negative. The true 
distances are the absolute values of these quantities. 
Since the next step is to equate these two distances, it is 
only necessary to show that these two juantities always have 
the same sign. There are two cases to consider, depending 
on whether the first component of the apparent position is 
positive (X>0), or negative (X<0). Letomrxt 0809) be the 
intersection of the X axis with the wavefront lane as it 
passes through the acoustic center. Taen (3.3) can be used 
to solve for X', namely 
M NI tc2''t CI) / 2.C1, 
Now consider the case where X>0. Then C1>0 also, and if 
leno) 2S positive, then 
a cT 
which implies that the wave arrives at the acoustic center 
before it arrives at the the x hydrophone. Therefore, since 
X>0, the plane at the acoustic center will intersect the X 
axis at a pcint beyond the x hydrophone, or X'>D, and hence 
Meet => (ET C2 + C3)7 2 C1 » 1 
AA EZ BES > 2 CT (Sines > 0) 
= C lt C i C3 » O 
=> (3.5) is positive. 
A parallel argument applies for the case of X>0, thus 


establishing that (3.5) and (3.6) always have the same Sign. 


Equation (3.7) is the result of eguating these two 
quantities. 
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For convenience, let K = 2V/D, and then apoly the same 
logic to the distances of the y, 2, and c aydrophones fron 
tne wavefront plane as it passes through tke acoustic 


center, to obtain the equations of (3.8). 


Ci ~= C. = C, - KT + K Ty. = 0 
(3.8) 

“Cj + e, _ C, - KT + K E. = 0 

“cy - C, + C4 - KT + K T, = 0 

e + e, + C, + KT — K e = 0 


The system (3.8) is four equations in four urknowns, and 
is the planar model's version of the equations (1.1). The 
unknowns in (3.8) are Cl, $S62pP 9G EGn TS However there is 
the additional constraint that Cl, C2, and C3 acze direction 
cosines, and therefore the direction cosines constraint 
(3.9) is a fifth equation, creating a a system of five equa- 


tions in four mi noq na, 


2 2 
1 2 3 (3.9) 


Generally there is no set of values (C1,C2,°3,T) which 
will satisfy ail five equations at once. This is because of 
the realities of a nonplanar wavefront and the presence of 
timing errors. The next section developes a method that 
produces set of values for the unknowns which is intended to 


satisfy those eguaticns reasonably well. 


D. MINIMIZATION OF SUM CF SQUARED ERRORS 


Let Ei, (1=1,2,3,4) be the value of the leit hand side 
of the i-th equation of (3:8): Then Ei measures the amount 
of error in the i-th equation caused by the chosen solution. 
It is impossible to have Ei=0 for all i=1,2,3,4. However 
soue ccmpromnise may be made. Specifically th2 compromise 
Chosen here is the classic minimization of (3.10), the sum 


of squared errors. 
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2 2 
Ej + E, + ES (3. 10) 


Minimization is to be done subject to tha direction 
cokinc#ecnstratat (3.9). Application of the Lagrange multi- 
plier technique [Ref. 5 p.55] calls for the minimization of 
tm Over all possible choices or Ci, C2, C3, T and 
lambda. 


È a2 S 2 
L= ALE - No E i) (391) 
i=l i=l 


Taking tne partial derivative of L with r2 spect to 7 
pelas (3.12). 


4 

OL 

A (m - -2K ( E. + E. + E. + E, ) 

OT à š IL 02, 


Il 


zs Las + ENS oum *] 


EUN: 12) EG zero ana Solving for T yizlás immedi- 
ately the appropriate estimate (3.13) of T, the ray transit 
time. 

T = E E I 
2 Į 2 3 4 (2159 

The partial derivative of L with respect to C1 is 

(510). 
or 
Cie 


È 


E EM iw ee 
1 2 3 4 all SH.) 


2 E + 2KT + K ( nil ) 3 


DINE eas used for T in (3.14), then (3.15) results. 





D _ ] 
iB - 2 Ca AD SEES 9-15) 
1 
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If the same procedure is used for the parcial deriva- 
tives of L with respect to each of Cl, TC? andy Coe cece 
are equated to zero, then the results are (3.16). 


P m AD c, | See WENN (3. 16) 


In order to solve for the Lagrange multiplier lambda, 


square both sides of the three equations of (3.16), and add 


-the resulting eguations together. Then use the sum of 
squares constraint (3.9) and solve for  lambia to yield 
(32175: 


3 
2 
A «NOM ele, ER 
l= 


Substitute (3.17) in (3.16) and simplify to obtain the 


appropriate estimate of Ci, namely (3.18). 


Les. 
1 


oes (3. 18) 


Ho. 
tl 





The choice of sign in (3.18) 1s determined by the fact 
that Ci is positive if and only if the sound wav? arrives at 
the i-th phone before it arrives at the c-phone, which in 
turn implies that (T4-Ti)>0. 


E. THE LEAST SQUARES METHOD 


In summary, the first alternative method for determining 
an apparent position starts with estimation of the ray 
transit time to the acoustic center by (3.19). Then the 
apparent position estimates are computed using (3.20). 


ENT (T. ST ee me (3: 19) 
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This method shall be hereafter referred to as the 'Ieast 
Does Method", or ‘L.S.' for short. The apparent advan- 
tages or the L.S. method are that: 

l. all four hydrophone times have egual veight in a 
Simple expression for the ray transit time T, rather 
than using an expression so heavily dependent on Tc 
as in the NAVY and NAVY_A methods; 

2. the differences of squared time values «hich appear 
in the soluticns of the NAVY methods are avoided in 
pue _ L.S. Method, thereby lessening tie tendency 
toward computational roundoff problems; 

3. the direction cosines already add to unity when 
Squared, requiring no arbitrary adjustmen: ; and 

4. cCccmputation of the initial angle allows cancellation 
of several terms, resulting in the simple expression 
(35.24). 


5 
A = arcsin ( > = T. ) RU de D _ È y? (3. 20 


F. BIAS IN THE LEAST SQUARES METHOD 


Unfortunately a potentially serious problem exists with 
puo n5. net hod. That concerns the consequeices of the 


assumption of a planar wavezront. The effect is difficult 
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to derive explicitly because it is difficult” oa ce na 
the true shape of a wavefront corresponding to a nonconstant 
velocity profile. However if it is assumed that a spherical 
assumption is more accurate the planar assumption, then the 
bias can pe estimated roughly, and then subtract ed from the 
original L.S. position estimate. This is SES a dierent 
problem because, as previously discussed, the four tasic 
spherical equations themselves have no exact solution. 
Nevertheless, as a rough estimate of the bias, the 
following procedure is used. First estimate be apparent 
position by the 1.S.  netkedé Then calculate the straight 
line distances from that position to each of th2 four array 
hydrophones. Divide those distances by the velocity of 
sound at the array to obtain the corresponding times. Use 
these times to recalculate the apparent position using the 
L.S. method ararn The difference between the >riginal and 
recalcualted L.S. positions roughly measures the error that 
would be made by the L.S. method when it is applied to the 
time values which correspond to a spherically spreading 
Sound wave whose source is in the vicinity of the original 
apparent position. Therefore this difference cai be used as 
a bias vector which can be subtracted from the original L.S. 
solution. This bias correction is the basis of the method 


set forth in the next section. 


G. THE LEAST SQUARES CORRECTED METHOD 


The secord alternative method for estimating an apparent 
position is as follows: 
1. calculate the arparent position P1 by using the L.S. 
nethcd; 
2. calculate the distances from that position to the 
four hydrophones, and convert them to times by 
dividing by the speed of sound at the array; 
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3. use the new times to calculate a new aparent rosi- 
EE USing the L.S. metnod ; 

4. calculate the difference vector P2-P1 (see figure 
3.1), and suttract it from the original position P1 
to ottain the corrected position P 

P = P1 - (P2-P1) = 2 Pi - P2 

5. finally adjust the cay transit time T calculated for 
the original position:P1, Dy Using the (proportional 
transformation: 

T = T OS R RI 

where R is the range to the new position ?, and R1 is 


the range to the original positior P1. 


Recalculated 
L.S. Estimate 
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Figure 3.1 Bias Adjustment for the L.S. Method. 


This method shall hereafter be referred to as the 'Lleast 
mares Corrected Method', ‘or L.5S.C. for short. The proper- 
mes of the LaS. C. method, like those of the NA/Y A method, 


are not well understood at this time. It is offered only as 
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an alternative which may combine the beneficial properties 
of the L.S. and NAVY methods, namely that it will: 

1. estimate a ray transit time value equally dependent 
on all four data times, thereby smoothing out exces- 
sive error in any one of the time values; and 

2. reflect the spherical wave assumption, believed to 
provide a more accurate description than the planar 
assumption. 

For targets at a range of 3000 feet, the l2ngth of the 
error vector P2-P1 varies from 0 to as much as 10 or 11 feet 
(see Table II). The error vector lengtns seem to be depen- 
dent cn koth azimuth and elevation angles of the target fron 
the array- These patterns indicate a potential for further 


investigation to relate estimation errors to suci variables. 


He | MAXIMUM LIKELIHOCD CALCULATIONS 


One shortcoming cf all methods described thus far is 
that none oí them can be used to produce an estimate of the 
underlying error in the time data values. The method set 
forth in this section is a first attempt to estimate that 
noise, and is again based on the assumption of a planar 
wavefront. 

Let Ti be the time recorded from by i-th hydrophone. 
Let Ji be the true time which, under absolutely error free 
conditions, would have been recorded at the i-th hydrophone. 
An assumption of Gaussian noise is now made, nam2ly that 

Td. ss Ud bn 
where the Ei are independent identically distributed normal 
random variables with mean zero and variance S2. Therefore 
the Ti (1215225375) are also normally distributed with the 


sane variance, but with means Di (i=1,2,3,4). 
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Now let Cj be the j-th direction COSTE A MN 
geometrically, using the assumption of a planar wavefront, 
Cj is given by (3.22), where V is the speed of sound at the 


array, and Dis the array dimension (30 feet). 


er 


j — er (3. 22) 


Letting U = U4, then (3.22) can be solved for each Uj in 
terms of U and Cj, as in (3:23): 


" y. f CD (EN 


The probability density of each time value Ti is given 
pv (25920) 


il 2C T, - E 
= —— —_—m 3. 24 
EE Jom s xP we ( ) 


Using (3.23) in (3.24), and multiplying the fourme neS 
ties together, the likelihood function (3.25) is formed. 


4 
-1 4 
e exp |, E Po - U + pc. /v)* \ 
PUE S i=] 


In (3.25) C4 has been used for notational :onvenience, 





and is defined to be zero. Then the log-likelihood function 


is formed by taking natural logs ofm(3:25),;, yedingmi 2c 


1 
SE 





L, = 2 In(2g7p wc 2 NS S 


4 
A | 
1 Lo (3. 26) 


Since the values of the Ci are to be direction cosines, 


the usual direction cosines constraint must be added to L1 


to form the Lagrangian function L2 given in (3.27). 
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3 
o I Lc - | i 
i-l -. 


The objective is to maximize L2 over tie possible 
ss es cn Cl, C2, €3, U, S and lambda. Mobbing S Lor the 
moment,  maximization of L2 can be acheived by minimization 


ere Given in (3.28). 











1 4 > E 2 
L = > c? > (T; = U + DE, /V) A ER A | 035029) 
1=1 1=1 
Wale partial derivatives of L to get (3.29) ind (3.30). 
OL 2D 
mo ea A = 2 Me. (3. 29) 
1 
4 3 
QL : | 
AS p m. - 4U + SS (350309 
Ou à 3 M TE 


sss - ID to zero and solve for U tolobtain (3.31), 


tke maximum likelihood estimate for the time at the c hydro- 


phone. 
4 5 3 
Safe 27; ye NN Si (3.31) 
3=1 J=1 
Eguate the three equations of (3.29) to zero, ani 


Ii each equatioóon by Ci respectively, to ob:ain (3.32). 


2 Dos 
== —a À—" AP — [= ` — 


Add the three equations of (3.32) together, and use the 


direction coSine constraint to obtain (3.33). 


(3. 33) 


Jai 


LA 
l! 
"Do 
H 
O 
I 
= 
Me 
O 
LA 
+ 
Jo 
Esc) 


Eguation (3.34) results when (3.29) is equat2d to zero. 
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T ee 
i 


E =o m i = Tore (3.38) 
V D 
$A -2 
Substitute- 33) into (3.34) to obtain (32.3 5)e eee 
Maximum likelihood estimate of the direction cosines, where 
U is the estinate calculated by (3.31). 


Tr WD 

C = A —— — — Í  Á==———= 
j 335 
i 3 SCH ( ) 
y e 


UJ 


As the reader will perhaps have noticed, tie equations 


OL (25:35) define each of the unknowns Ci interms of all 
three unknowns. Such a structure suggests tiat eee 
be used as an iteration function. That is, 1: reasonatle 


initial values are used for the three unkrowns in the right 
hand side of (3. 35) then new values are produc2d. Repeat 
the process using the new values until the answez stabilizes 
within acceptable tolerances. Aithough convergence to the 
correct solution is not guaranteed, the method has never 
failed for the equations of (3.35). Unlike the L.S. method, 
(3.35) does not have any known closed form solution. 
Returning to the standard deviation S, take the partial 
derivative of (3.27) with respect toS to Jec MP 


7 | 
F Dand M eee B D 2 
_ S 3 Lu Ure (3. 36) 
Multiply (3.36) by S3 and solve for Sè? to get the 
maximum likelihood estimate of the variance, given by 
(BRS) 
4 3 
2 i i D 
s = FA aa (3. 37) 
1 | i=1 
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I. THE MAXIMUM LIKELIHOOD PLANAR METHOD 


In summary, the third metnod for estimation of an 
apparent position is as follows: 

1. let U= T4 initially; 

Zi. use the L.S. method to calculate the initial values 
Mord. gsimngvt5e3i2])* 

See doe Uesand Ci (1-21,2,3) in the right hand side of 
(3.35) to obtain new estimates for Ci c 

ü. recalculate JU, using i 31). 

5. reiterate steps 3 and 4 until tae values Ci (i=1,2,3) 
and U converge within acceptable toleranc:2s; 

calculate S=, Using (3.37); 

7. calculate the estimated apparent position relative to 
the c-phone, using (3.38); 


X = V U = = 
I c! EN V UC, ^ V U C3 (3, 38) 


8. lastly translate this solution and its c)rresponding 
time estimate to a solution and time relative to the 
acoustic center, using (3.39), where R and Rc are as 
deramed Dy (4): in Chapter ti. 


X = X + D/2 YSP yY +D c cs 
š i (3. 39) 


R 
T UR / R. 


This method shall be herearter referred to as the 
"Maximum Likelihood Planar Method', or M.L.P. Tor Short. 
Originally the hopes for this method were ratner high, espe- 
cially since it was the first method to produc» a variance 
estimate. However subsequent experience with the method 
Ica tes that It probably suffers significantly from at 
least two factors: 

1, the planar wavefront assumption probably builds ina 


position bias as in the case of the L.S. nethod; and 
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2. the variance estimate is inflated since part of the 
noise being measured is due to the inadequacy of the 


planar assumption. 


J. COMPUTATIONAL EXANPLE 


For a quick comparison, the three methods leveloped in 
this chapter are now applied to the example which was used 
in Chapter: The EXACT results caiculated previously are 
included in Table III for comparison. 





a CO G T T WW Q" 

| TABLE III | 

Single Example Comparison of | 

Planar Wavefront Methods | 

Transit Time Elev. Angle | 

Method I (secs) A (degs) | 
task 0.67529319 15.22798 

ES 0.67527149 15 326 72 | 

Me cere 0.67529007 15. 10437 | 

EXACT 0.67527043 15.27 002 | 

Apparent Position Estimate | 

Las. ( 1003.809 , 3019.751 , 866.126 ) | 

| eS ace ( 1006.089 , 3018.266 , 868.235 ) i 

| M. 3. P. ( 997.722 , 3023.685 , 859.355 ) | 

| EXACT ( 1005.966 , 3017.899 , 868.555 ) | 


- 


ee E HR A E eg, i ee CC N de ee TM. RÍA 


As can Le seen easily, the M.L.P. method fares rather 
‘poorly in all regards, even worse than the L.S. method. 
This will be confirmed by the evaluations made in Chapter V. 
Also worthy of note is the apparent tendency of the L.S.C. 
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NAO Lo Correct the L.S. method back toward the exact 
solution. The evaluations of Chapter V will confira that the 
L.S.C. method almost always yields a better solution than 
Aor ginal L.S. solution. 
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IV. A SPHERICAL WAVERFRONT MODEL 


A. THE SPHERICAL WAVEFRONT ASSUMPTION 


All the models developed in Chapter Iii ire limited 
primarily by the assumption that the wavefront is planar 
upon arrival at the hydrophone array. In this charter a 
model is developed under the assumption that the wavefront 
is spherical upon arrival Mat the array: If the sound 
velocity prcfile were constant with depth then the spherical 
model would be exact. This 1s of course not the case, but 
it is suspected that the wavefront is better modelled asa 
sphere than as a plane because that small piece of the wave- 
front which passes through the 30 foot cube spanned by the 
array is locally spherical. That is because every part of 
that piece travelled through approximately the same regions 
of water, experiencing the same general raybendiig patterns. 

The spherical assumption is accurate if and only if the 
speed of sound is ccnstant over the ray path, and conseg- 
uently the originai four spherical equations apply once 
again. They were: 


2 2 | 2 


t9 


( X - D/2 ) + (YDI) IECUR ADU = 
^ _ x (4.1) 
( X + D/2 ) + ( Y = D/2 e + (2 + D/2 E = ve 
(x + 0/2 )2 + ( Y +. apup 02020) a 
( X * D/2 1 + ( Y + D/2 17 + ( 2+ D/2 » = vir 


It has been stressed previously that there is no exact 
solution (X,Y,Z) satisfying all four equations (4.1). That 
is because the time values on the right hand side correspond 
to the reality of a variable velocity profile. However, if 


the spherical wavefrcnt assumption is to be accuzate, then a 
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constant velocity is the assumed case, and any innacuracies 
in the time values are regarded as due to timing errors 
Oni, Therefore, under the spherical assumption, if the 
true time values Ui (i=1,2,3,4) were known and substituted 
into (4.1), an exact solution to the overdetermined system 
would be realized. In that case a solution to any three of 
the equations would also be equal to that unique exact solu- 


Elon, im particular, the NAVY solution of Chapter II would 





be the true solution. Therefore in terms of the coordinate 
System referenced to the c hydrophone, X would be given by 
(4.2). 
2 
D V E 2 
However, X is also given by (4.3), wher» C1 is the 


HN Citron Cosine along the X axis of the vector from the c 
hydrophone to the sound source. 


au Mm | (4.3) 


The time value of interest at the moment is J4, the time 
at the c hydrophone. Therefore let U = U4 for clarity of 
presentation, and then eqguate the expressions in (4.2) and 
(4.3) in order to solve for U. If this same l»gic is also 
applied to the similar expressions for the distances to the 


y and z hydrophones, the results are (4.4). 


2 2 
U, = Yu - (2DUC,/V) * (D/V) zu 2.3 (4.4) 


The expression (4.4) will be useful in the development 


of the model of this chapter. 


B. LEAST SQUARES MODELS 


A direct approach might be to apply the l2ast squared 


error technique to the spherical equations, in a manner 


4g 


paralleling that used on the four planar equations of the 
L-S. modelin i Chapter ma However the fcrmulae and eçua- 
tions that result are exceedinyly complex, involve fourth 
degree powers of the data values Ti, and have thus far 
defied all solution attempts. Therefore this id:a was aban- 
doned in favor of the maximum likelihood approach which 
follows. | 


C. MAXIMUS LIKELIHOOD COMPUTATIONS 


As in the M.L.P. model, Gaussian noise is assumed for 
the time data values. Therefore 
Tl = Ui + El 178 1722255 
where the Ei are independent identically distributed normal 
random variables with mean zero and variance  S?. The 


density of each Ti is therefore (4.5). 








' S 1 = 2 
GC = >° EE ( Ti = U, ) | i= I US (4.5) 


Io form the likelihood T func tron, multiply the four 


densities together, to obtain (4.6). 


E 3 | 
L = (—— exp = > (malo : + (7 -U) 


` 5 





Forn the log likelihood function EE EE 
rithms of (4.6) to obtain LI of DE 





3 
i B 2 z 2 
L ln(27[) - 4 1n(S) - A Leon + ED | (4.7) 


In order to maximize L1 with respect to C1, C2, C3 and 
U, it is sufficient to minimize L2 of (4.8), where a substi- 


tution for Ui has been made using (4.4). 
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? 
de 


3 
2 L 
= E n di ( OS ` 
L, ( 4 U) > | 5 U (200 m + (D/V) 


gre 


| (4.8) 


cosines constraint and forn 
Bueeagrangian function L of (4.9). 


3 
es AE d (4.9) 
i=l 


For notational ccnvenience, let Ki be given by (4.10). 


Now add the usual direction 





ie 2 
en Ü = (2DU C; / V) * (D/V) i fgg (4. 10) 
Take the partial derivative of L with respect to Ci to 
get (4.11). 
ƏL , D (T; - Ki) E : E 
E -2 x. V E AC; mo 10,3 : 


Simpliry (5.11) and equate to zero to yield (4.12). 


cau i = 1,2,3 
i —— 7 M M 
V K, 
oE 


Multiply the three eguation in (4.12) 





(4. 12) 


DIE Tal, 273) 


respectively, and add then together. Than use the 


constraint on the sum of squares of the direction cosines to 
obtain (4.13). 





3 

mc 

8 DU i Ni 
i=] VA 

(4. 12) to get the max. mum likeli- 

hood estimate for the direction cosines, as in (4.14). 


Substitute (4.13) into 
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Lee ee 
n (4. 14) 
r = 


3 
TEES 
f isa 
- J = 
j=1 Sj 


Now take the partial derivative of L with rzspect to U, 
as in (4.15). 


3 Es. 
OL = > da Sa 
JAS ). SF eu -bewavy|  — ONSE CE 





Equate (5.15) to zero and solve for IU to obtain the 
Maximum likelihood estimate of the time value (4.16) at the 
c hydrophone. 

3 E 
= = 28 des e (4. 16) 
j=1 We 


SC, K 
AS 

l -= 
b 


Finally take the partial derivative of L1 in (4.7) with 





respect to S. Equate it to zero and solve to get (4.17), 


the maximum likelihood estimate of the variance. 


3 
2 : 2 "EM - 
S D [ie a ) Í (T, U) | (05017) 


D. SOLUTION BY A MODIFICATION OF NEWTON'S METHOD 


The solutions given by equations (4.14), (4.16), and 
(4. 17) once again form a set of equations which would seen 
to be solvable by naturai iteration as in the M.,.P. method. 
Unfortunately this time the technigue fails to converge. A 
mathematical tool is needed which is stronger chan natural 


iteration. What is used is a modified four dimensional 
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Veeoton OL WNewton's Method {Ref. 5 p.47] to seirch for the 
roots of a set of four equations. 

The cbjective is to determine the values C1, C2, C3 and 
U which satisfy (4.14) and (4.16). For further notational 
convenience define the values Mani Nas in (4.18) and 
(4.19). 


Së = EL (4. 18) 


yas) di 


Given those definitions of M and N, then tie equations 


whose roots are desired can be simplified to (4.20). 


T. -«4JK. = 
B , i v i (4.29) 
C: n; (C4, C, C. ,U) = — i= 
K. M 
Vi 
U = h : 
4 04:05:04, 0) 
1 — N 


Now define error functions as in (4.21). "base evaluate 
the amcunt of error in each of the equations (4.20) for any 
porwof walues for (C1,C2;,C3,U)- 


g. = C. - h. j DR 
L là , 
AS (4.21) 
Let G be the four dimensional coli mn vector 


(ons 2eq3 5:94) *.. Finally let GP be the matrix of partial 
derivatives of G, as in (4.22). 

Newton's Method in four dimensions says that if X(n) is 
a four dimensional vector holding the current approximate 
Fr s TCO C 3 and U, then X(r+1) will be an improved 


answer, where X(n+1) is given by (4.23). 


DIS 

















SED CEN (m 27) 
ERR oc. 
Cer NOME 
GP (C. Co C. U) = OC, OC, 
Inn o5 
Dr o. 
S. s 
oc, Oc, 
X = X - [ce axo ] ; G (x) 
—+1 — =n (423) 


Of course it is necessary to calculate the derivatives 
held in the matrix GP in order to use (4.23). Tiose deriva- 
tives are given in Appendix B. 

Unfortunately when multidimensional versions of Newton's 
Method are applied, there is often a tendency for the method 
to converge slowly, or even diverge. This is because it 
tends to overshoot the best answer for each iteration. To 
alleviate this problen, a modification is made to the 
method. At the end of each Newton iteration, PELOE to 
proceeding with the next iteration, a Golden Section Search 
[Ref. 6] is performed to find the best possibl2 answer in 
the direction of the new iterative solution. Specifically 
the line in four dimensional space from X(n-1) of the 
previous iteration to X(n) of the present iz:eration is 
searched for the best answer. The definition of the 'best' 
answer is that point along the search line whith minimizes 


the sum cf squared error functions, namely (4.24). 


Y | 
g. š TENA 
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mencu rent iterative solution X(n) is thei given the 
value orf the minimizing point resulting from the Gclden 
Section Search. Then the next iteration of Newcon's Method 
is performed, along with another Golden Section Search to 


find the next iterative solution X(n*1). 


E. THE MAXIMUM LIKELIHOOD METHOD 


In summary, the fourth and final alternatire method to 
estimate apparent positions is as follows: 

1. let U= T4 initially; 

2. use the L.S. method (&e321) to initially estimate the 
values ofkeD (jal,2, 3); 

RENE 0): 

4. initialize the Newton Iteration counter : I - 1| ; 

5. calculate the values Ki, HM and N in accordance with 
equations (48.10), (5.18) and (5.19); 

6. calculate the error function vector  G(X(I)), using 
(4.20) and (4.21); 

7. calculate the derivative matrix GP(X(I)), using the 
results in Appendix 5; 

Sue pinvert GP (X(1I)); 

9. calculate the new Newton estimate X (I+1) from (4.23); 

10. perform a Gclden Section Search along the line 
between X(I) and X(I+1) to find the point which mini- 
mizes (4.24); 

11. let (+1) be egual to the minimizing point found in 
the previous step; 

12. increase the Newton iteration counter : I = I+1 ; 

13. reiterate steps 5 through 12 until the values 

ES S 276980) 

Converge within acceptable tolerances; 

14. calcualte the estimate of S? using (4. 17).. 


DO 


This method shall hereafter be referred to as the 
‘Maximum Likelihood Spherical Metro d OLI O 
As will re seen from the evaluations nade in Chadter V, the 
M.L.S. method 1S apparently the only alternative to consis- 
tently rival the performance of the currently used NAVY A 
method. It has the additional advantage that it estimates 


the error present in the time data values. 


F. A COMPUTATIONAL EXAMPLE 


This latest method, M.L.S., ls now applied to the same 
example considered in Chapters II and III. “Or EISE 
comparison, Table IV lists the results of using all the 
methods. The error vector lengths are the distarces of each 
position estimate from the EXACT apparent position. 

Since this is crly one example, this table is not 
presented for the purpose of any broad conclusiois. However 
it is of interest to note that in this example 

1. the M.L.S. method outperforms all others, including 
the NAVY_A method; and 

2. in many ways the original NAVY method out»erforms the 
adjusted NAVY_ A method. 


G. VARIANCE ESTIMATICN 


In the computaticnal example, the time data values used 
were exact since the methods of Appendix A coild be used 
with the known linear velocity profile. Therefore the 
appropriate value for variance in the time valies would be 
Zero. For this error free example, the estimates shown in 
Table V were obtained for the standard deviations of timing 
noise, using the two maximum likelihood methods. 

In the case of this one example, the spherical assump- 
tion is apparently an improvement over the planar, since the 


M.L.S. estimate o£ error is only 4% of the M.L.P.. estimate. 
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E A E 
| TABLE IV | 
| Single Example Comparison of All Methoús f 
| 
| Method Esas a A Veet or 
NAVY 0675 208/59 15.26356 0.4128 | 
| NAVY_A 0. 67321027 15.26461 0.4840 | 
| tens 06732999 1522795 j 733 | 
| Tris, 0. 67527149 15.26372 09207 | 
| il. L.P. 0:216 1521990 7 15. 104 37 13.64) | 
EL. S. 0.67527C49 15.26677 0.4117 | 
| EXACT 0. 67527043 1527227902 | 
| Apparent Position zstimate | 
NAVY (10057957 3017.874 SG Senso) 
| NAVY_A ( 1006.087 30 187259 Sod. 29589 | 
iS. (100 609 3019:7773 866.126 ) | 
| SC. ( 1006.089 20127266 Sedo | 
UE (9941 122 3029056895 959023595) | 
| USL. S (1006. 195 BO le. 190 865: 1753) | 
EXACT I ( 1005. 966 30177399 5685555” | 
E" 1 ] iut. | 
That is regarded as an improvement because.the correct 


answer is zero. The nigher M.L.P. estimate is indicative of 


the inflation due to the planar assumption. 


The tracking range which provided data for this study 
records time values to seven decimal places. Tierefore the 
Standard deviation estimated by the M.L.S. m2thod is of 


particular interest, Since it indicates errors in the 


seventh decimal place even when there is no ercor present. 


Since the data was actually error free in this exampie, the 


sui 


— aaa 


IABLE V 


Maximum Likelihood Error Estimates 
for the Example Problen 


Ven (defteg Ke 


Hodel Method Est. Std. Deviition 
Planar MLP: 5.517 EID secs. 
Spherical MOSS 221917 Ez secs. 


es ee Ce ee Ob ———————— ee M ee ee ` mmm 


estimate is a measure of the variation induced by the spher- 
ical wavefront assumption. A broader discussion of error 


estimation will be undertaken in Chapter VI. 
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V. EVALUATICNS OF MODELS 


A. GENERAL 


Mere ate many sources of errors in the ov2rall hydro- 
phonic tracking problem. These include, among others: 

1. errors in estimation of the sound velocity profile; 

2. inhomogeneity of the velocity profile over time and 
horizontal displacement; and 

3. possible errors in measuring the positiois, and the 
angles of tilt and rotation for the hydrophonic 
arrays. 

This study focuses on those errors which occur during 
the computations preceding the ray tracing procedure. TO 
evaluate the performance of the methods developed in 
Chapters II, III and IV, it is necessary to control strictly 
the ray tracing procedure. Only in that way can the differ- 
ences found between methods be attributed to the differences 
between models, and not to any source of error oi1tside those 
methods. 

It was originally hoped that the various methods might 
be compared by applying them to real tracking data. However 
it was found that the overali tracking problem aad too many 
large sources of error to allow the methods to demonstrate 
any differences. Therefore the methods were con pared under 


a more tightly contrclled simulated environment. 


B. SIMULATION SCENARIO 


Two different simulations are used to compare the six 
different methods. Both use a basic scenario similar to 
that of the computational example explored in Chapters II, 
III and IV. That example assumed that: 
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1. the velocity versus depth relationship is linear, and 
given exactly Ly: 

V = 4840.7 + 9.03314 ° DEPTI ; 

2. the acoustic center of each array is at a derta of 
1300 feet; 

3. the hydrophone arrays are ail level, and their X, Y 
and Z arms are parallel to the respective coordinate 
axes of the tracking range. 

Under these circumstances the methods se: forth in 
Appendix A can be used to compute the exact values for the 
hydrophone times, ray transit time and elevation angle fora 
sound wave emanating from a source at any specified loca- 
tion. Therefore when the methods are appiied to those exact 
times, the resulting estimated positions can be comparea to 


the known true position. 


C. SIMULATED ERROR FOR TIMING DATA 


The models developed in this study were designed to 
improve position estimation, especially in the presence of 
errors in the timing data. Lacking any better model at this 
time for timing errors, the simulated environment includes 
an assumption of Gaussian errors for the hydrophone times. 
Therefore realistic timing data can be simulateì by adding 
to each exact time value a random guantity of normally 
distributed error. The mean of the error is assumed to be 
Zero. The variance was estimated from real tracking data, 
using the variance estimating property of the M.L.S. model. 
The data from one tracking run was used, involving six 
hydrophone arrays and 733 position estimates (se? Chapter VI 
for data selection details). Each position from the 
tracking run produces one estimate of the variance. The 
variance value chosen for use in the simulations was the 
median of the 733 variance estimates produced by the M.L.S. 


method. That value was 
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Sé = 9. 1204 b= 12 secs? . 
That is the same aS a Standarad deviation of 
S = SE EE secs . 

Each of the two types of Simulations was ru: four sepa- 
rate times. Each run was done with a different specified 
variation. Those four distinct error conditions are delin- 
eated in Table VI. 


——_ —————T_—— 


TABLE VI 


Simulation Error Levels 


RUN LEVEL VARIANCE STD. DEVIAĽ ION 
1 Zero 0.0 2.0 
2 low S 5-149 93.402 E-T 
3 medium 9.12 5-212 s. 02 B-E 
4 high 9.12 E-10 EEO2 DS 


[ee  À MÀ to Gët E deefe Dec o, ia 


TE SEW SCH 
| 
| 





D. SINGLE ARRAY SIMULATION 


In the first simulation, the intent was to compare the 
methods pairwise, so aS to determine which me:hod is nore 
likely to produce the more accurate estimate of a given 
sound source position. One thousand positions were chosen 
at random. Each position was 3000 feet from the array. The 
positions were uniformly spread over the surface of a sphere 
of radius 3000 feet, centered at the array, truncated above 
by the water surface (depth 0) and below by the depth of the 
array (1300 feet).  ZIhe methods set forth in Appendix C were 
used to assure that the random positions selected were 
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uniformly distributed over the surface area of the truncated 
hemisphere. 

Each of the 1000 randomly chosen source positions vas 
then processed as follows: 

1. calculate the exact hydrophone times, using the 

methods of Appendix A; 

2. add to each of the four exact times a random value of 
error at the specified level (zero, low, medium or 
high); 

3. apply each of the six methods to the hydrophone 
times, generating six different apparent dositions; 

4. apply the ray tracing procedure to each of the six 
positions, using layers that are 25 feet thick, and 
utilizing the known linear velocity profile, thereby 
producing six different estimates of the sound source 
cca tion. 

5. compare the six different estimates pairwise to see 
which method in each pair produced the estimate 
closest to the true sound source location. 

The comparison being made is that one method is consid- 
ered preferrable to the other if it more frequertly produces 
the more accurate estinate. 

The layer thickness of 25 feet was selected order to 
Simulate the actual procedure at the tracking range which 
provided data for this study. However, as discusséd in 
Chapter I, when the velocity profile is smooth and known 
exactly, the process is very robust with res»ect to the 
thickness used. The thickness values 1, 10, and 25 were all 
tried, with virtually no changes in any of the comparisons 


between methods. 
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E. SINGLE ARRAY SIMULATION RESULTS 


Tables Vil and VIII contain the results of the single 
array simulation. Each tabled entry represents :he fraction 
of time that the method of that row produced a better esti- 
mate than the method of that column. For example, in Table 
VII, the L.S. method outperformed the M.L.S.~ method in only 
5.8% of the 1000 trials with low error values. 

In a one tailed test that one method is detter than 
another, these binomial proportions are significant at the 
0.05 level if they exceed 0.526. Symmetrically, one method 


is significantly worse than another at the 0.05 level if the 


proporticn is less than 9.474. For the 0.01 level the 
corresponding critical values are 05:53. and 0. 463 
respectively. 


The results indicate that: 

1. as previously clained, the NAVY_A method usuaily 
outperforms the unadjusted NAVY method; of particular 
interest is the case of zero error whith actually; 
compares the relative ability of each method to 
produce the exact answer when given the exact times; 
in those cases the NAVY_A method does extremely weli 
against the NAVY method; 

cera error Conditions the most successful 
performer is the M.L.S. method, since it always has a 
favorable (greater than 0.5) comparison fraction 
against ali other methods; the M.L.S. fractions vary 
little over the four error levels; 

3. under all error conditions, the worst performer is 
always the M.L.P. method; 

4. spherical methods consistently outperform planar 
methods; and 

5. increased error levels tend to lessen the distinction 


between methods; in the zero error case comparisons 
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TABLE VII 


Single Array Simulation Results 
for Lower Error levels 


uU 
| | 
x 
| ERROR LEVEL : ZERO | 
| NAVY NAVY A L.S. SA MIS Masa i 
| NAVY mom .950 26905 ¿987 . 420 | 
NAVY A .989 2950. 4328 389 Wau | 
LaS: |. 050 - 050 JO | 987 |. 057 | 

Le Bee 2324 222112 |. 950 (987 . 386 | 
M.L.P. 2013 . 013 201.3 z043 Bnp | 

| MLS. .580 . 566 .943 .65Uu |. 987 | 
| | 
ERROR LEVEL : LOW | 

NAVY NAVY A L.S. Leo nea Mo LaPa bd. Las, | 

NAVY . 165 |. 952 .648 . 987 . 390 | 
NAVY A .835 |. 952 |. 693 2987 .405 | 
Todd S 2045 . 048 .048 .987 .058 | 
SE 6352 22202 4,952 a VB} . 369 | 

Ma. La. P. 01% . 013 .018 . 013 «043 | 

We L A 604 « 595 « 942 a631 . 987 | 


| 


are in the 


interval 


(0.01,0.99), 


error case that interval 
(0. 3770 16 3)32 


F. DOUBLE ARRAY SIMULATION 


In the second 
compare the 
method is more likely to 


closely in the two array crossover problen. 


Simulation, 


methods pairwise, 
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the 


this tine 


while in 


intent again 


the high 


was 


1S narrowed considerably to 


to 


determining which 


produce positions which agree more 


This is not the 


, 





| TABLE VIII 
| EU mi Simulation Resuits 
| gher Error Levels | 
l ERKOR ISVEL : MEDIUN | 
| NAVY NAVY A L.S. DISCO MALOS. 1 
NAVY . 464 +019 59 29597 4H 35 | 
| NAVY_A .546 .821 .566 .987 .437 | 
i LS 2:151 . 173 e451 A . 177 | 
| tee ce 44841 - 434 2919 NS - 450 | 
| TSP. 013 . 013 -029 eo V3 0.3 
Mo. 0.565 21503 219253 25510 . 987 | 
x PFRAQOR evel + “HIGH 
NAVY NAVI A L.S. Laso. PS AOS 
| NAVY . 486 + 544 2439 25027 . 431 | 
NAVY A .514 . 546 .516 .559 . 431 | 
| I. e456 . 454 . 457 252m . 400 | 
fp eo. Coe 0 1 . 484 «54:3 2502 2027 | 
Ber». .438 . 441 - 443 . 438 . 373 | 
IUE SS. 569 . 569 .600 2573 2927 | 
| 


same question as that addressed by the single array simula- 


tion. The two estimates produced by any one m2thod may be 


very close to each other, and yet be far away from the true 
position. 

For the double array scenario two arrays are used, 
rated by 7500 feet, koth at depths of 1300 feet., 


are parallel to the 


sepa- 
The arms 
of each array corresponding coordinate 
system axes. Once again 1000 positions were raniomly gener- 


ated in a uniform manner, this time over a 3 dimensional box 
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Figure 5.1 Double Array Simulation Configuration. 


running crossways between the two arrays. The box is 1300 
feet deep, 5000 feet long and 1000 feet across (see figure 
5.1) 
Each of the 1000 randomly chosen sound source positions 
in the two array Simulation were processed as follows: 
1. calculate the exact hydrophone times for the first 
array using the methods set forth in Appeidix A; 
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2. add to the exact times some random error at the spec- 
ified level; 

DE Repeat steps 1 and 2 for the second array, using a 
new set of random error values; 

4. apply each of the six methods to the time values from 
both arrays, Producing Six differen: pairs of 
apparent positions; 

5. apply the ray tracing procedure to bo:h apparent 
positions in each of the six pairs, utilizing thie 
known linear velocity profile, prođucing six pairs of 
estimated sound source positions; 

TOR Sach of the Six palrs of positions, cilculate the 
distance between the two positions in the pair; 

7. make pairwise comparisons of the six different 
distances, to see which method in each pair exhibits 
the closest agreement between its two position esti- 
mates. 

This time the ccmparison being made is that one method 
is considered -better than another if the positions it 
produces agree more closely more often than those of the 
Other method. 


G. DOUBLE ARRAY SIMULATION RESULTS 


Tables IX and X contain the results of the double array 
simulaticn. Each  tabled entry represents the fraction of 
time that the method of that row produced a pair of esti- 
mates which were in closer agreement than the estimates 
produced by the method of that column. Significance 
criteria for these proportions are the same as for the 
Single array simulation. 

These results are similar to those of the single arrav 


simulation, because they indicate that: 
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| 

o TABLE TZ | 
Double Array Simulation Results | 

for Lower Error Levels | 

ERROP LEVEL : ZERO | 

NAVY NAVY_A L.5. et e oie ©. OS | 

NA VY . 366 . 9069 29 75 .968 2212 | 
NAVY A  .634 . . 735 -968 2212 | 

l LaS: .011 2011 <O 29939 -014 | 
LISSCITUVZ7U « 2606 2989 - 968 . 195 | 
Dee 052 -445 2032 033 | 
Niele $785 . 788 . 936 . 807 2967 | 
ERROR Deven, -SELON x 

NAVY NAVY A L.5. L. S S£ C S “eb. 2. O | 

NAVY . 448 . 989 ~0 07 2 9912 . 30 | 
NAVY_A 0552 2989 -684 + 952 2312 | 
eo .011 oon -011 . 560 |. 018 | 
L. S C mn s 2 16 TO 29 T24 | 
M.L. p. WZ 00 8 - 048 .440 .048 |. 036 

| Me LoS: J690 -688 -986 .746 .962 | 
Me — d PPP — 


the NAVY A method outperforms the original unadqdju 
NAVY method, although the difference is not signifi- 
cant in the higher error level cases; 

the most successful performer is consistently the 
MeL: S: method; 

spherical methods almost always outperform planar 
methods; 

increased error levels tend to lessen the distinction 


between methods. 
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TABLE X | 
Pee te | 
ERROR LEVERS: MEDIUM | 
NAVY NAVY A L.S. L.S.Co Bebe Ee re Leone 
. 498 . 841 «513 .808 . 443 
2902 ~ 842 «242 .809 - 440 
Sale? = oe 2 159 «525 . 160 
.487 . 458 .841 .808 . 433 
4 MET . 191 2475 212 . 144 
-957 . 560 . 840 «2957 . 856 


ERRORELEVEL : HIGH 


NAVY NAVY_A L.S, pi ue. PESE 

. 487 . 556 .535 .479 - 442 
.513 .558 .499 .481 .440 
444 . 442 .445 BER .422 
.465 . 501 .555 .480 . 438 
.521 .519 2549 .520 | 428 
.558 . 560 .578 ET .572 


AAA A o op QA BEE ve no DA cess A 


| 
| 


are however some 1lndications from these results 


different from those of the single array simula- 


tion, such as: 


1. the worst performer was consistently the L.S. metnod, 
rather than tke M.L. P. method; 
2. in the high error case the M.L.P. method is equal to, 


Or even marginally better than, every other method 
except the M.L.S. nethod; 
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3. the performance of the M.L.S. method 15 noticeariy 
better under error free conditions. 

The last two inferences are perhaps the most inter- 
esting. The maximum likelihood approach was intended to 
estimate and account for Saussian errors in the timing data 
values. Hence it is really not very surprislig that the 
M.L.S. method does well with error prone data. BES 
interesting to note that the M.L.P. method also does well 
under high error levels, even though it probably suffers 
from a bias due to the planar assumption. 

Of even greater interest 1s that tne M.L.S. method seems 
to Le at its best when compared to other methods under error 
free conditions. This result was unexpected, and indicates 
that the M.L.S. method not only handles timing errors well, 
as was intended, but apparently aiso does an eve1 better job 
of approximating the elusive exact solution to the original 
four spherical eguations of (4.1) and (2.1) when the exact 


time values are available. 


He LIMITATIONS ON INTERPRETATION OF RESULTS 


The results of both simulations seem to im31y that the 
M.L.S. method outperforms the currently used NAVY A method 
on any randomly chosen sound source position, with or 
without timing errors. These are encouraging results. 
Nevertheless the reader is cautioned that thes» tests are 
just simulations, and like all simulations, must make 
assumptions which cannot fully reflect the reali: y of actual 
hydrophonic tracking conditions. The most important assunp- 
tions made for these simulations are: 

1. the sound velccity profile is known exactly, and is 
linear; 

2. the errors in the timing values from any hydrophone 
are normally distributed with mean zero, and are 


independent of the noise in any other hydcophone; and 
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3. the parameters of the simulation were iixed; Fe, 
example the single array simulation used a fixed 
range oí 3000 feet, and both Simulations used arrays 
at 1300 feet of depth, with fixed orienta: ions to the 
range coordinate system. 

The first assumption is probably of little -onseqguence. 
It not only greatly facilitates conputations, but also helps 
to isolate the initial angle and time estimation € problem 
from the unrelated errors involved in the velocity profile 
estimaticn procedure. 

The second assumption is somewhat more troublesome. 
Errors may not be Gaussian at all, or if they are, the mean 
ray not be zero. Unfortunately each positioi estimation 
involved only four equations, and therefore did not allow 
for estimation of more than four paramenters. Tierefore the 
error mean, being a fifth parameter, could not be estimated. 
Also the errors of any one hydrophone may very well not be 
independent of the errors of the other three phones on that 
array. Fortunately these concerns are offset somewhat by 
the results of the double array simulation, wherein it was 
found that the M.L.S. method was at its best when there was 
uomnorse at all. 

Also the error type and level may depeni on other 
factors, such as the target's range, elevation and azimuth 
angles from the array. This highlights the concerns of the 
third assumption. There is considerable room here for 
future work concerning the dependency of results on such 
coupiacating factors. 

Lastly it should be pointed out that the simulations 
make comparisons only on the binomial basis of better versus 
worse in 1000 trials. The magnitudes of the actual differ- 
ences are iynored. It is possible, though perha»s unlikely, 
that while one method marginally outperforms a second method 


in most trials, in all the remaining trials the :irst method 
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is much worse than the second. There is also room here for 


further work. 
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VI. CONCLUSIONS AND RECOMMENDATIONS 


A. ESTIMATION OF TIMING DATA VARIATION 


One of the primary purposes of this study sas to esti- 
mate the amount of variability in the timing data being 
Wecorded during actual tracking runs. This problem was 
addressed by the M.L.P. did -P-5 elaxinum Jlixelihood 
models. The M.L.P. model, as previously discuss2d, suffered 
from a cias due to the planar assumption, was the poorest 
estimator of positions among all the models, and produced an 
inflated variance estimate. Therefore the spherical model 
M.L.S. is used to estimate the data variance. 

The variability that is measured ty the M.L.S. model is 
Made up of three componerts. First there is the variance 
induced Ly the spherical assumption. Then there is the 
variance caused by the seven decimal accuracy used when 
recording the data. Finally there are the errors inherent 
in the physical process, due to such factors as hydrcrhone 
variability or malfunction, local distortions of the sound 
wave, and inexactness of the water column which estimates 
the speed of sound profile. The last two sourzes of error 
together make up the variability that is involved in the 
time values which are ultimately used in position estina- 
tion, and is therefore the variation that is to be 
estimated. 

If it is assumed that the variability indiced by the 
spherical assumption is independent of the data variability, 
then the M.L.S. variance estimate is the sum of those two 


variances, or 


2 = 2 + 2 
O ais O sph O ti me 
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Therefore the data variance can be estimated Ey first 
estimating the variability induced by the spherical assump- 
tion, and then subtracting it from the M.L.S. estimate of 
the variation. 

The M.L.S. estimate” Was obtained by applying tre. .l-o- 
methcd to the data from a tracking run at the Nanocse 
torpedo tracking range on May 6, 1980, That run invclved 
position estimation Ey several different hydrophone arrays. 
The run made several thousand position estimates, 733 of 
which were at depths of 100 feet or more and involved 
targets not more than 4700 feet from the sensing array. The 
depth liritation was imposed to avoid the excessive corpli- 
cations caused by the radical changes in the velocity 
profile above that depth (see figure 2.1). The maximum 
range linitation imitates the data validation procedure at 
the tracking range, where positions farther than 4700 feet 
from the array are discarded. 

The tracking data yielded the following range of esti- 
mates for the standard deviation of the data noise, using 
the MEL.S. model 


MAXIMUM VAIUE 2409 5-5" (secs. 
0.95 QUANTILE 1.18 E~5 secs. 
MEDIAN VALUE 3.02 E-6 secs. 
0.05 QUANTILE 4.11 E-7 secs. 
MINIMUM VALUE 3. 13 E-8 secs. 


For an overall estimate of the noise, the median value 

was used, so that: 
2 = 307 5E 
Ü nis 

The variability induced by the sphericai assumption was 
estimated Dy applying ene sine s procedure to perfectly 
noiseless data in the idealized environment of the single 
array simulation of Chapter V. This was done for targets at 
ranges of 1500 to 450C feet, at 500 feet increments, with 
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1000 randomely chosen targets at each range. The results 


are collected in Table XI. 


A ee EE O M E E AZ NA AA s stiess E e EE bf i ee m e 


TABLE XI 


Inflation of Error Estimates Induced 
by the Spherical Mode 


Standard Deviation Estimates from Exact Tim» Values 
with Known Linear Velocity Profile 


RANGE MINIMUM  Q(.05) MEDIAN Q(.95) MACIMUM 
1500  2.72g-10 2.84E-8 2.08E-7 3.13E-7 3.29E-7 
2000 2.122-9 5.88E-8 2.26E-7 3.145-7 3.30E-7 
2500 3.90E-8 1.09E-7 2.34E-7 3.14E-7 3.31E-7 
3000  8.35z-8 1.50E-7 2.37E-7 3.15E-7 3.24E-7 
3500  1.21E-7 1.59E-7 2.37E-7 3.13E-7 3.24F-7 
4000  1.47E-7 1.65E-7 2.39E-7 3.16E-7 3.25E-7 
4500 1.46E-7 1.69E-7 2.42E-7 3.14E-7 3.25E-7 


E ia EECH 


hum BA IEEE 


| 
| 
| 


Table XI shows that the median and maximun inflation 
values are reasonably independent of target range. The 
minimum values vary somewhat, but only for range values 
below 3000 feet. This represents a very stable situation 
overall. Therefore the inflation due to the spherical 
assumption is estimated by the median value at a range of 
3000 feet, namely: 

2 = (fee) ene e p. 
Ü sph 
Combining these two estimates, the varian>e estimated 


for the timing data is 
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2 = 2 _ 2 
O time Onis O sph 


( 3 02 EDO” Toa > K Ee ET 
(359 dI me e 


As can be seen, the error induced by the spherical model 
is less than  J0* ot thes error estimate. Therefore 
when it is accounted for by subtraction from the  M.L.S. 
variation estimate, the final variance estimite changes 
little. 

The estimated value indicates a standard error in the 
6th decimal place. With a typical speed of sound value of 
4880 feet per second, this represents a position differen- 
tial Tof about 

4880 e 3.011E-6 = 05015 Leora. 
This estimate is guite iow, indicating that the time values 
being recorded are sufficiently accurate. 

There is considerable opportunity for addi: ional work 
determining the relationship, if any, between the time vari- 
ance and other factors such as angles of elevation and 
azimuth of the target from the array. 

There is also the problem that the time 7ariation is 
likely to be array dependent. For example, consider figure 
6.1 , wherein the standard error estimates from the actual 
tracking run are plctted versus the range of the target. 
The plot does not indicate that there is any simple rela- 
tionship between range and error level. Howev2r the plot 
does show a bunching pattern. | When the error estimates are 
plotted separately for each array, then tne bunching pattern 
becomes clearly associated with the individual arrays. 
Consider figures 6.2 and 6.3, where the separate plots have 
been made for four different arrays. It is still noggmeleass 
from these plots whether the principle effect is due to the 
individual arrays, or the ranges of the targets. However 
some level of array dependency seems likely, indicating a 


need for additional investigation. 
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DATA FROM ALL ARRAYS 
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Figure 6.1 Error Estimation Yersus Range of Target. 


ee ‘CHOICE OF METHOD 


Clearly all indications are that the planir wavefront 
models, L.S. and M.L.F. are not candidates for use as posi- 
tion estimators. Furthermore the hybrid model L.5.C. is an 
interesting improvement, but never really performs well 
enough ccmpared to the M.L.S. and NAVY models. 

The original NAVY model is usually outperformed by the 
adjusted NAVY A  nmethcd. However the differences are not 
always significant. 

The spherical model M.L.S., on the other hard,  consis- 
tently outperformed all other methods during the simulated 
evaluations. It wculd seem that M.L.S. is tie model of 
choice. It loes the best job of handling normally distrib- 


uted errors in the data. But that is not tae strongest 
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Figure 6.3 
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argument for its use. A more important, and SurpBe songs 
argument in its favor is that when the exact, error free 
times are used, the M.L.S. apparent position estimate will 
usually produce the most accurate estimate of the true 
apparent position. This is the desired overall result, so 
that the sensitive ray tracing procedure will be affected 
Minimally ky apparent position estimation errors.. 

There are nevertheless several notes of caution which 
should be considered before embracing the M.L.S. method 
wholeheartedly. The first caution has been stressed before, 
namely that these conclusions were arrived at under the 
idealized conditions of the simulations scenarios. The 
second caution, also previously stressed, is tha: the actual 
magnitudes of the differences between position estimates has 
been ignored. It is conceivable that while one method 
always produces a better estimate than another, the differ- 
ence between any two fosition estimates is acceptably smail. 

Lastly there is the caution that the M.L.S. model 
involves a complicated iterative procedure which uses 
considerable computer time. It is probably too slow a 
procedure for use with ‘real time't analysis during the 
execution of tracking runs. For real time tracking the 
NAVY_A method currently in use is probably preferrable due 
to its simple computations. 

However, for post run analysis, and also >ossibly for 
calitration of the hydroph one arrays, the M.L.S method is 
recommended as being a more exact and more robist position 


estimator than those methods currently in use. 


C. RECOMMENDATIONS FOR FUTURE INQUIRY 


Several recommendations have already been made for work 
needed to estimate the effect of suitable independent vari- 


ables on both timing errors and the bias in certain metnods. 
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In addition there exist at least two other areas for 
possibly fruitful investigation. 

The first area concerns the interplay betwzen methods. 
Specifically, the binomial comparisons of Chapter V show 
that even the worst methods are better than each of the 
other methods at least part of the time. Hence it is 
possible that the best method overall would be a suitable 
combination of methods, wherein each method is used where it 
is most effective. For example, even though the M.L.S. 
method has been indicated as the best method for any 
randomly selected position, it may be consisteitly outper- 
formed by another method under certain circumstances, such 
as extremely high or low elevation angles. 

The second area for possible work addresses the question 
of how to next improve upon the existing models. It is 
herein suggested that the next improvement in modelling 
would be a method which is based on a linear velocity 
assumpticn. As figure 2.1 demonstrates, a linear velocity 
profile is a reasonable approximation for most depths. This 
would be the next logical step above the constant velocity 
assumption which is associated with the spherical nodels. 
Most of the mathematical basis for such a model in contained 
in Appendix A. Possibly a suitable set of equations could 
be developed involving the hydrophone times and reflecting 
the linear velocity assumption. If so, the least squares or 


maximum likelihood techniques might provide useful results. 


81 


LINEAR VELOCITY PROFILE THEORY 


All computations in this study are rade under the 
assumption that the sound velocity is directly related to 
depth'in a linear manner, and is known exactly. Under these 
circumstances many closed form results, not othe: wise avail- 
able, can ke obtained and used in those computations. 


Suppose that the velocity profile is given by 


y= VO + V1 z 
where VO and V1 are known constants, and Z is the depth 
variable, measured down from the water's surface. In this 


case it is known [Ref. 4] that the path of a solnd ray from 
a signal source to a hydrophone will have the shape of an 
arcot a- cinsslei The center of that circle will be some- 
where above the surface of the water. The vertical place- 
ment of that center is determined by tne value of z at which 
the speed oí sound equals zero (see figure A.1). Although 
that depth is negative and is not really a depth at all, it 
nevertheless has geometric meaning. 

Consider the vertical plane containing ec cie 
center, the sound source and the hydrophone. set h be the 
variable which measures the horizontal position in that 
piane. Let (h,z)=(ai,a2) be the position of the hydrophone, 
and let (h,z)=(pi1,p2) be the position of the sound source. 
Then C2 given by (A.1) is the Z coordinate of the circle 


center. 


C = 
2 
ET (A. 1) 


what must be found is Ci, the hl coordinate oi the circle 





center, and r, the radius of the circle. To solve for these 
values, the eguation of the circle is used, evaliated at the 
the two known locaticns, yielding (A.2). 
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Figure A.1 Circular Ray Paths for a Linear Velocity Profile. 
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The left hand sides of the two equations of (A.2) can be 


eyuated, and solved for C1, leading to (A.3). 


_ (P4 * a) | (P, = a.) SH m 
S 2 Die ED». (A 
i 1 
When this value of C1 is substituted into the first 


equation of (A.2), then the result is given by (A.U). 
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The circular arclength between the hydropk»ne and the 
sound source is easily computed. Unfortunateiy the velocity 
of sound does not stay constant along that arc. Therefore 
in order to determine the amount of time (T) “equired for 
the sound ray to travel the ray path, the effect of the 
velocity must be integrated along the arc. This is dcne by 
(A.5), where S is the arclength between the two points, and 


V(s) is the speed of sound as a function of post tion on the 


arc 
S 
x 8 dz 
" VIS (A.5) 
O 

eee aie) the sound source position corresponds to an 
arclength of s = 0, and the hydrophone position correponds 
to arclength S = S. In [ Ref. 4] this inteural is shown to 


be equivalent to (A.6). 


SH 
n 28 E EE) 
Vi cos ia) 
20 


In (A.6) A0 is the angle of elevation of the ray path at 
the sound source, and A1 is the elevation aigle at the 
hydrophone. The antiderivative in (A. 7) can be used to 


solve (A.6), leading to the ray path transit time expression 


in (ARS): 
da a pl LAA) 
cos (a) Gos fa (A.7) 
cos (A. ) l + sin(aA 
r = — n ndi — (2.8) 
1 coc IE: sin(A})/ 
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If the elevation angle at any point along the arc is 
denoted ty A, then (4.9) relates the angle to th: derivative 
of Z with respect to h along the ray path. 


dZ 
tan (A) = SIR (A .9) 





Gn acit differentiation or the equation of the circle 
yields (A.10) as another expression for the same derivative. 
T h = Ci 
dh B Z- C (A. 10) 





Equating these two derivatives leads to (1-11) which 


relates the elevation angle and the position (h,z) on the 


arc. 
h - C1 
C e o (A. 11) 
2 
Pron (A. 11) a Simple geometric argument 2oroduces the 
cana tions 15 (A. 12). 
= C. h - C1 
cos (A) = r sin (A) = FRE (A. 12) 


First let A, h and z be equal to (41,al,a2) SCENE E WÉI 
and then let them equal ({A0,p1,p2). Then substitute both 
expressions into (A. 8). The result is (A.13), the desired 
expression for the ray path transit time in terms of the 
positicns of the sound source and the hydrophone. 


a. = C r + p. - C 
me =. in Ge E s (A. 13] 
1 P = C, LN a; - C4 
To Summarize, if a ray path is to terminate it the three 
dimensional position (X1,Y1,Z1), and the sound source is at 
ENE LEA perform tne following steps in order to 


calculate the exact ray path time and elevatio1 angle that 


would correspond to a linear velocity profile: 
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23 
ere 


use (A.14) to translate the tnree dimensional fosi- 
tions to the two dimensional coordinates of the 
previously described vertical plane, with the origin 
at the water surface directly above the ena of the 


ray path; 





(A. 14) 


calculate C2, Cl and © uSing {A. 1), Arlen s 
calculate the exact ray path transit time T using 
(A. 15); and 

calculate the exact elevation angle A at the hydro- 


phone, using (A.15). 


a = m 
2 2 
A = ——_ —— —. —T 
arccos ( ) (A. 15) 
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APPENDIX B 


PARTIAL DERIVATIVE FORMULAE FOR NEWTON'S METHOD 


The central tool used by Newton's Metnod in the develop- 
ment of the M.L.S. model in Chapter IV is the matrix GP. It 
is the matrix of first order derivatives of the error 
expressicns gl, g2, g3 and g4. Those derivat. ves are set 
forth in this appendix. 

mek = mcl,CZ2,C35,U), then the error functions are 
defined by (B.1), 








Sh i Sh 
g; (X) > SAA i S (B.1) 
P M 
ui 
T, - D M / Ç 
SHX) = U - 
1 - N 
where the values Ki, M and N are the functiois given by 
MZ) . 
2 
E - D 
K, = U ( el ) + 5 (5.2) 


Z 
H 
UJ E PIG 
uo 
d 
E 
LJ 
d 
C. 


Recall that GP is the matrix shown ir (8.3). 
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GP (C1, C5, C4, U) - Ci Oc, E 
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Then given the definitions above, the folloring meras 


tives are the result of straight forward, though tedious, 





(B.3) 


differentiation, and are offerred without detailed proof. 
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(B.4) 


(B5) 


(moo) 


(B.7) 


ECR AULA a 

The first three diagonal elements of GP are the deriva- 
tives of each gi with respect to Ci (i=1,2,3) aid are jiven 
EE EC, 


OM 
39; x= WE - VK,(T; - JK; ra Pe 


- 





EOIMNUDLA 2 
The cff diagonal elements in the first three rows and 
columns of GP are the derivatives of each gi with respect to 


qnd are given by (B.9). 
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ECRMULA 3 
The derivative of each gi with respect to U is given by 
(0:10). 


E S = OM 
NON IT, Ay com (B. 10) 


OU 
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O` 
G= 
< 
< 
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R 
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x 
P: 


FORMULA 4 

The first three elements of the last row of GP are the 
derivatives of g4 with respect to each Ci, and ire given by 
ESTI): 








On (vr, - DM) E O" ay) 
99, OS OC; (B. 11) 
Oc; v (1 - n) È 
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FORMULA 5 
The last row, last column of SP is tho ri oo n K D, 


with respect to U, and is given by (B.12). 





ON EN (Om 
ER 3 


V (19 ND) 
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APPENDIX C 
UNIFORM SAMPLES ON A TRUNCATED HEMISPHERE 


The single array simulation oí Chapter V required a 
random sample of positions in space uniformly distributed 
over the surface of a truncated hemispnere. Th» hemisphere 
is tc be of radius r (3000 feet) about the acoustic center 
na nydrophonic array- The truncation of the overall 
Sphere 1s due to the fact that the upper portion of the 
hemisphere is above the water surface, and the lower half is 


telow the sea botton. 





w = Y COS(E) dH 


dud 
\ 
| a ay AS AAA MI ciere, A A amne qui) 
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Figure C.1 Hemispherical Geometry. 


Let E be the variable denoting angles of eleration above 
the horizontal. Let H be the variable denoting horizontal 
azimuth angles about the center of the hemisphere. Consider 


a small piece of hemispherical surface area bounded by the 
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elevation angles ET and (z1*dz), and by the azimuth angles 
Hi and (H1+dH) (see figure C.1). If Y is the borr ortal 
width of that piece, then W is given by (C.1). 


w = r dH = E COS LE) on (C.1) 


If Al is the area of that piece, then Al is i pproximated 
by HEHE 


A, = wdE =  r cos(E) dH dE (C2) 


Now suppose that ( 0 «^ EBEN AS Emix ) where 
Emax iS the elevaticn angle of the top of the truncated 
hemishpere. Then the ratio (A1/A2) of two different areas 


at elevation angles Ei and E2 is given by (C.3). 





A, m cos (E.) dH dE cos (E, ) 
eS h P (C3) 
A, E cos(E,) qi d cos (E) 
If n1 and n2 are to be the (relative) Sampl> sizes from 


the two areas Ai and A2 respectively, then uniformity of the 


sample requires that (C.4) hold. 


A ONUS 


(C4) 
The combination of (C.4) and (C.3) implies (2.5). 
cos(E_) 
n = n — _ 2 
2 J. cos (E, ) (C.5) 
Letting El = 0, then the relative sample size n2 for 


area A2 is given by (C.6), where n is the relative sample 


Size at the base of the hemisphere. 


B5 = n cos (E,) (96) 


Now the differential probability of drawiig a random 


position that has elevation angle E is given by (C.7), where 
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Nis the total sample size to be drawn from the surface on 
the sector bounded by the azimuth angles H1 aná (H1+dE). 


: mcn: 
fo (E) poc S y os (E) (Com 


Therefore, i£ K ls defined to be the constan: ratio n/N, 
then the corresponding cumulative distribution is given by 
(Ge. 8) o 


E 


= á = Kem E) “~ 
UE 1 COS e) de (C .8) 


Bormcusmnigpestactevatione angle, Emax, th» cumulative 
distribution function must equal one, so that (C.9) holds. 


D nE (C.9) 


As a result, the constant K is determined by (C.10). 


di 


K = —MáM— 
SU y) (C. 10) 


Therefore the complete cumulative distribution function 
is given by (C.11). 


u sin (E) 
Fp (E) 7 Sin(E ) (C. 11) 
max 


Now the inverse frobability transform can be used. If J 
is a uniform random variable on the interval (0,1), then let 
Ee TA, 


E = arcsin Cr 
sin(E  ) (C. 12) 
max 


Now choose an azimuth angle H randomely aid uniformly 


over the interval (0,270). Then (X,Y,2) given by 


r cS5s(H) cos (E) 


r sin(H) cos(E) 


r sin(E) 


9s 


is the corresponding position in spherical coordinates. 
That position is a random position drawn from a population 
or positions uniformly distributed across the sirface ofa 


truncated hemisphere with maximum elevation angle Emax. 
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SINGLE ARRAY SIMULATICN COMPUTER PROGRSA 
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FIR3VAHH ha tts RANA s o *RNOO0rrAana NH 
(ye OI. COM OP D Do DeAH n bri)» > e ~A 
ZUHA m H QUO HAUMAMNOE HAM YEAH se fu a O z 
Tora Oz A> HRAAHN BAMUMHN A DI e> (t tq OA E O a 
AH BOO nar uc EH Cn SA HH x y (n ° : 
MNFOUH A QE ZU MH ZH < «OD mtt NA c =a 
c= Ei = Die OO -«amP»IOOum ZO bic DI O "MO 
OMHNIYWH D HGHH Sm mHOt q MC Em A nea Uu) AW HD 
m= “mU Q BP Ati NOH HHAax O AGOO0OVOO HA OHP GO HO 
WHEGHO WN TO Da bis QAH mm fa sHHHHHHK az r Om ei Hu) 
QO n BHHUNMH DH MAND MAaAMnuununn pa AZA PRA NH 
Oa m Ey ett A NODO He ' HIHHHHHH E SN OA 
eërbabitäe O MOOG e UO QUO zu) HAL) +3 O0» OQO3 nO 
ed tH O (GEI sbb) Het OA sw Dip ex Azf €v—Owe at 
EA tF E= A Emu) z ORG QO == mm meme Ho ° **¿—OMO WA 
m. O E H MAH Goobit/O m=kH- c m= 43E-i D OMG HO DOMA < DA mn 
O m A AM Ot AGW xXiotrinindHeduU) NH GH HA OO “N00 « (hg 
hs eun A ei X2 Y HHDAO0Hz0 S Hh ANDAMAN Ba CO III LA I N mu 
ZH na HnOoxAa WHO UAHS BS 02239 HA enm emn edu) 
na tu m MAOZ HAG DO U) F3 FJ ea ea m e ea ea c XiE4E4O II 
HN GHA bis eme QUMAMUNHR HX HSD9DDIONII HA IN M P3 PR FA EUH HAG 
ti TH OH-XAHA GO DT mini TH BHOOCOOOO mz AOMGHH mo 
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AAO 


ANNAN 


AN 


100 


0000 
i undi 
O O N 
è ee SD 


an) 


200 


"p = 4840.7D0 
VV 2 = 3314.D-5 
V 2 VV(1) +VV (2) *A (3) 
SET UPPER LIMIT ON ELEVATION ANGLE, IF NEC? SSARY 
IF Í DE cm a (3) JO TOO 
HIMAX = 1.2566D0 
COTO 20 
PHIMAX = DARSIN(A(3) /R) 
GENERATE RANDOM VALUES FOR ANGLES OF 
ELEVATION AND AZIMUTH. 
FK = 1.D0/DSIN (PHIMAX 
CAB SCHEER 
CALL GGUBS (SEED,M,REL 
CONVERT ANGLES TO 3-D COORDINATES 
DO 11 I = 1,M 
EL(I) = REL (I 
AZ (I) = RAZ (1 
PHI = DARSI 4 L 1) FR) 
P(I,1) = R*DCOS(PH1)* cos MN : 
P(1,2) = R*ICOS(PHI)*DSIN(AZ (1) *P1*2.D0) 
P (1,3) = A(3)-R*DS IN(PHI) 
CONTINUE 


COMPUTE EXACI HYDROPAONE CPES UNDER TETE NEER 
VELOCITY PROFILE ASSUMPTION. 


CALL TCOMP (P,AC,M,VV, TC) 
CALL TCOMP (P,A,M,VV, T) 


GENERATE AND ADD ERRORS TO TIMES. 


NNOISE = 4*M 
ac GGNMNDWSEED ,NNOIS Eos TU Pr) 


DO 111 TE RM 
12 


nO 122 M = Wee 
NOISE - STUFF (E) 
NOISE = NOISE*S DEV 
TN (I,J) = T(I,J)+NOISE 
K = K + 1 
CONTINUE 
CONTINUE 
FORMAT (2X,'* ( SAMPLE SIZE =  ',D,' )') 
RITE (6,35 1 M 
WRITE (6,999 
THO OUTER LOOPS RUN THROUGH ALL PAIRS OF TIE SIX 
POSSIELZ POSITICN ESTIMATING METHODS. 
DO 1111 METH1 = 1,5 
KMETH = METH1 # 1 
DO 1222 METH2 = KMETH,6 
IF (METH1 .EQ. META2) GO TO 1222 
FORMAT (2X 'METHCDS USED ARE: ') 
WRITE (6, 200) 
METH = MPTH1 
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G TI O [IDEPOUTTHESEIOPITIMPLEHENT TEE METHODS 


Ir (METH 1) GO TO 131 

ALL Posi y (TN, vi MP rT O TES 
FORMAT NAV ONCORRE TED") 
il 387) 


IF 
TEST) 

NAVY, CORRECT: D') 
DE 
EST) 
S, JNCORRECTED') 


HU 


TE 


mH 


LA tu ` Ca Lech Ce E ` sch Oz CH ` sech De CIE dad Oe 


IF 5) GO TO 135 
P E UA, PES 
T4 I 


I TEST UAR) 
. LIKELIHOOD, PLANAR') 


em~ HI'O 


MG 


HH Hime 
rÜ < 
mi 


1 
) 
FIH1,METH2 
T FIND METHODS *,14,* AND/OR *,I4) 


IH 
HEA OSAMA SIONS IO. SIAM SIO, 


ION ESTIMATE 
e AND MUST BE 
GE COORDINATE 


8:3 : 30 
ST (I,3 


ARA Qs 
Aint VOI 
RESTI Spi ii 


trj prn Ej 
Mette Ae OPMOWOreowoYreZOwmOoOrzezowmworzoworm 


O el 

O O 

Gi 

H JU = 

H 

= 

G 
tori 
e Ed Ed 
UN un 
— HH 

NHAN 


CCNVERT APPARENT RUE RER Re EK BY 
RAY TRACING PROCEDURE. 


Ceo CETTPESI,PILTESEA,M,VV) 


COR UTATE DIPEERENCE BZIWEEN THE 15T POSITION 
ASIN AT AND THE TRUE POSITION. 


IF 2 METH .NE. METH1 ) GO TO 160 
O DIE I) = 0-DO 
DIFP LS) = “DABS (TEST (I) -TC (1,4) 
it 


156 J = 1 
DIF(1) = DIF(1)+(P7 (1,3)-P (1,9) )*» 
CONTINUE 
CONTINUE 
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OQ GOO GO O 


160 


167 


165 


UI LYLY 
OOO 
NWO 


OVO IN 
OOQ N 


— fN 


METH = MERN 
GO TO 130 


CULATE DIFFERENCE Be PEEN THE ZA eo eee) 
MATE AND THE TRUE EEN ANDAS 
HEIST POSITIONED P ENENGE; 


I 
© 


[33 CJ [3 UO HA zi) F+ O 
UN Hmm Ga 


I — 
- 
a 


1,3 
) = DIF(1)-(P7(1,3)-P(1,3))**2 


FT(I) - DABS (TEST (1) TC (1,4 
HM Ti 0. DO y “Go td 165” ?) 


es DO ) GO TO 166 


H HOA VHamu HH 
HHO ODA 


HA WHG I 


ES 
CT, 
s + 
I (NYZS) /DFLOAT 
ACTION OF TIME 
IN P 
IN T 


zi 
tri 
H 


tozz PHA BHH 
tt HH 


tia OO sas 
ZKROO OO why 
owm zz HHH 


SIS 


DOUBLE ARRAY SIMULATION COMPUTER PROGRAM 


WY = 
(D PS e P4 
=O = E 
HAGUR 


NDOM SAMP. 


RÁ 
NS 
TC 
TO 
d 
E 


Eno «atm 
D Ofd'ta cz 
FQ O4 C E41 x EA 


Z1 FONNI 
eG) Qm B eH m 

A NH 
HRW m 
6 O > E S4 ed 


Xi o3 ex ea t- n3 
CAU) C) iG ex CD ed 


QOOOOQOQOQOQOQOQQOQQOQOOQUOQOOQOQOOOQOOOOOUO VU O 


V 


DEO SCUND PROFILE 


O PER 
DiS GIVEN BY 
Cho: 


IH 
AN 


H> MAS 


THE POS 


ERRORS ARE 
BUTION IS NO 


UND WAVE ARRI 
ION SDEV. 


ARE COMPUTED 


4 


EA 
SC 
M 

CM 
RT 
AT 


AOO fu E E 
EIU. 20 vw 
THUMM 1 H0 0 
SW) ALBIN 
=> m EC O 
H ORE Mam 
N YELU ZAU 
ANE Lo = 
AAT MATOZ 
OOXDODODH 


A 110 EM EH 
mun =, DO O O 
"igi oss otaztn 
HAZ Boo 
NA eh m4, H 
HHEAN D> uU 


TAMAOHS 


SEH POSSIDLE Pi TIR OF 


EE T K METHIMERWII]  HETH O NYES ,NNOLIS E;NARRAY 


DIMENSION RP(1000),STUFF (4000),NAM1(6),NAM2(6) 


` `. ` & 








DATA NAM1 
DATA NAM1 
DATA NAM2 
DATA NAM2 


OOUO00 


DEPTH 


RANGE, 


DUESDIPALIZE VALUES FOR SAMPLE SIZ&, 
AND ERROR STANDARD DEVIATION. 


Elio) 


DEV. = 


SE SID; 
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ere 


C 
3 


AQQ 


AQQ 


AAAN 


— N) 


AN 


CL 


Wn 
UNO) 


0 


EU 


A 
-RANDON ERRO 
x 
( 


tx 


T UP ARRAY POSITIONS, AND SOUND VELOC I I EE b sp: 


<< NNN a — sech 
NA UNA UN) 


le n 


( SAMPLE SIZE = #71 500) 


(SEED, M , RP) 
I) 

(XYZ- 0.5D0) *BOX (J) 
CONTINUE 
CALL GGIB 


BO 33.1 = 
XYZ 


SR 
CONTINUE 


CONPUTE EX RIVAL TIMES, ANDES 


R 
FS. 


CALL ICOM 
Coco. Sol 


Ero 


CONTE 
CONTINUE 
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125 


131 


132 


133 


134 


MOANA (O 


144 


OQ O00 


SET UP LOOPS TC RUN THROUGH ALL PAIRS OF M2 THODS. 
DO 1111 METH1 = 1,5 
KMETH = METHI, t da 
DO 1222 METH2 
IF (METHI . EQ. "oz Gor Te). 1222 
FORMAT(2X,'METHCDS COMPARED ARE 1. 2e 
FORMAT(2X," AND 2. *,44,2X,AU 
WRITE (6, 10 urn). Naz (A METH 2 
WRITE(6,20) NAMT(METH2) , NAM2(METH2 
WRITE (6,999) 
METH = METH 
NARRAY = 1 
SET UP ARRAY BEING USED 
A y = Al 1) 
A(2) 2 A112 
2(3) = A41(3) 
DO CE-I = 13 
E IJ) Dn J) 
CONTANÓ E i 
CCNTINUE 


Sree sus ROULTNES TO PERFORM ESTIMATION METAODS. 


IF (METH .NE. 1) co TO 131 
ALL POSNAV (T,V, M, PEST, TEST) 
GO T9 150 
IF (METH .NE. 2) GO TO 132 
CALL POSN VC (T,V,4,PEST, TEST) 


GO TO 
(HETH NE. 3) GO TO 133 
ALL POSLS (T,V,4, PEST,TEST) 
GO T2 150 
IF (METH .NE. 4) Së rO 134 
LLL POSLSC IT, ,PEST,TEST) 
GO TO 150 
IF (METH .NE. 5) GO TO 135 
AIL POSMLP (T,V,M,PEST,TESI,VAR) 
GO TO 150 
IF (ME NE. 6) ES TO 136 
CALL POLS (T; AS PSA E) 


WRIT 2 (6 137) METRI ME TH2 
FORMA Li ‘Land PIND METHODS ',I4," AND/OR ',I4) 
GO TG 

Tr (uA KAYU FO. 2 k GO TO 152 

T POSITION ESTIMATES ARE IN LOCAL À 


2RA 
ATES, AND MUST BE TRANSLATED TO TRACKI 
YSTEM COORDINATES. 

I 


M 

T(I,1)'- BEST (I, EIS HU 

T 1.2) = PEST(I 2) + A1 (2 
- A1 


LE 


a 
NG 


BEST 
CONTINU 


CORRECT APPARENT POSITIONS BY RAY TRACING. 
CAP R OCET PESI PI, TESTA, M, VY) 


101 


OQ GO OO 


152 


145 


OGOGO 0) 


OOO 
amb — 
Wn 
UIO 


OQ000 


160 


167 


-N 


OOO = wech 
OOD0O N) 
OIDO =N 


GO ON TO SECONL ARRAY 


NARRAY = 2 
n = A2(1 
A(2) = A212 
A E = A2(3 
DO 77 I= 1,M 
Se Am ba J) 
a E , 
conti ngi 
CONTINUE 
GO TO 130 
TRANSLATE 2ND ARRAY APPARENT POSITIONS TO 2 ANGE 
COORDINATE SYSTEM, AND THEN RAY TRACE. 
DO 145 I = 1,M 
PEST(I,1) = PEST (I) + a2 (1 
PEST 172) = FEST(I,2) + A2(2 
PEST(I, 3) = A2(3) - PEST(I,3 
CONTINU 


CALL TRACE (PEST, P2,TEST,A,M,VV) 
COMPUTE DIFFERENCE VECTORS FOR 1ST METHOD 
IF ( METH .NE. METH1 ) GO TO 160 
DO 1 = 1,M 
DIF “DO 
ge ) "StF (4) (21 (1,3)=22 (150112 
= + è 
CON j i i 
CONTIN 
GO ON TO SECOND METHOD 


META MET A2 


N 

Ej U de UN 
Ud AC ll j 
lO.» 


I 
I) 
56 
Ir 
IN 
UE 


ARRAY = 1 

GO TO 125 
CCMPUTE DIFFERENCE VECTORS FOR 2ND METHOD AND 
COMPARE TO DIFFERENCES FOR 1ST METHOD. 
NYES =0 
DO 166 I = 1, 

E MI D "ix I)-(PTQI,J)-P2 (00 655 

CONTINUE 7 (JI i 

IF í DIF (TI) .GT. 0.D0 ) GO TO 166 

YES = NYES + 1 

CONTINUE 
WRITE (6,999) 
FRACT = (DF QAT(NYES) /DFLOAT (M) ) 
FORMAT(2X,' FRACTION OF TIME METHOD #1 ') 
FORMAT(2X," IS BETTER IS £ A 
WRITE (6, 300 
WRITE (6,310) FRACT 
WRITE (6,999 
WRITE (6,998 
WRITE (6, 999) 
CONTINUE 
CONTINUE 
FORMAT (ax "+. o A 
FORMAT(2X;" di 
STOP 
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APPENDIX P 


COMPUTER SUBROUTINES FOR TIME CALCULATION AWD RAYTRACING 


OQOOOO001000000 00 


AAD 


22 
11 


OQA AAA 


O 


WE 
We 


SS SI Oo O O O O O O O O OI 9 990949 3 9D 


SUERNOUIINE TCOPENP)A,0,VV,1I) 


THIS FORTRAN SUBROUTINE COMPUTES THE EXACT TIME 
REQUIRED FOR A SOUND WAVE TO TRAVEL FROM ICS 
SOURCE (VECT OR P) THE FOUR HYDROPHONES ON AN 
ARRAY WHOSE ACOUSTIC CPE Sp EE EE E EE 
VECTOR A. 
THE METHODOLOGY USED IS SET FORTH IN APPENDIX A 
OF THE THESIS. THE BASIC ASSUMPTION IS ONE OF 
à LINEAR VELOCITY PROFILE, WHOSEZ COEFFICIENTS ARE 
GIVEN BY THE VECTOR VV. 
INTEGER M,I,J 
DOUBLE PRECISICN P(10 3) 4A (3) g VE 42) 7(1900,4) 
DOUBLE PRECISION AA(4 í ei) n 
SET UP COORDINATES OF FOUR HYDROPHONES ON THE ARRAY 
po A 

Gr Se S » 15. DO + AI 

= = : E 
HIE. ioo Am 
z 2 + 

CONSENTE 
AA(i1,1) = 15.D0 + at’ 
AA 2°29 = 5 DO + A 

373) = -15.D0 + A(3) 


CADCULATE QUANTITIES C1, C2, R, AND THE TI! ES T(*,ü) 
C2 = -1.D0*(VV (1) /VV (2) ) 

LOCP THROUGH THE M SOURCE LOCATIONS 

DO 33 I= 1 


e 


) 
LOOP T ROUGH | EE e EE Es CE SOUPCE 


SS 51 = SBT ( (2 (1 gi} 78a 1)) **2 
AA (J 2) ) **2) 
C1 = seed, ieee: Gs J,3) #*2 
+2.DO*C2*(AA(J, 3)-P2) 
C1 = C1/⁄(2. D0%B 1) 
La P On ngo pu 
Ea E RE 61:84) RCI) 
T(I,J) -» (DLOG(T(I,J)))/VV (2) 
CONTINUE 
CONTINUE 
FETURN 


END 
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SUBROUTINE TRACE (P, PT,T,A,M,VV) 


THIS FOETRAN SUEROUTINE TAKES AN APPARE 
ESTIMATE P AND CONVERTS TO A N ACTUAL 
IMATE PI BY RAY TRACING. THE RAY TR 
NEAR VELOCLIY PRO EIL WITH CORECTIE 
OR VV. OTHER INPUTS ARE THE ARRAY L 
OR A AND THE RAY TRANSIT IFE IMM 
ER, T. “THE RAY IRACING LAYERS 200 
THE FIRST Jud LAYER BEING THE 
IMMEDIATELY ABOVE THE ARRAY. 


NZH 
E Lia CG 
i ty (3 
O Lo 
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0200] 
pH AG O 


bd H 


EI 1000 0, 37 và (3) 
ih j 


ty 9 tU ta 


DEL 25 DIU 
LOOP THROUGH THE M APPARENT PO TTIONS 


ñ 


ESA TO VALUES POR 
Do 2) 
bs t6 au lt 
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= 
= 
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to, 
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NU H P 
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N+He 
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GO TO 20 
DEPTH 


H 
Noo 

“ty 
cl 


Da sch Fäi MINN D> 
i ro 
Hiro 
= 
Pilly TANU Pra Ril? ùs +€ NOt 
KO 


Í 
«4 


«T ty (7 


203 0 UDOOOUO HO < COON? il 
~ e t 


MWAH es Le Le Le lé 
O be 


O 
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O 
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ADJUST FOR OVERSHOOTING IN LAST LAYER 


= VDT- TKI 
= y*(D7-1(1)) 


ADJ APPARENT POSIT TONNO GER GC I IA POS IT ION 


PEIRONA 


DS 
DH 
UST 
PT (1,1 
PT (I, 2 
PT - 
T INU 
URN 
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COMPUTER SUBROUTINES FOR POSITION ESTIMATION METHODS 


2? 
Ti 


SS SS SIS EE EE EE ES EE 
SUBROUTINE POSNAV (T,V,M,P,NT) 
INTEGER MW, I 
THIS FORTRAN SUBROUTINE IMPLEMENTS THE ORI3 INAL 
UNADJUSTED NAVY APPARENT POSITION ESTIMATION METHOD. 
INPUT ARE THE HYDROPHONE TIMES T FOR M SOUYD SOURCE 
POSITIONS, AND VELOCITY V AT THE ARRAY. 
OUT2UT ARE M APPARENT POSITION ESTIMATES P ALONG WITH 
THE CORRESPONDING M RAY TRANSIT TIMES NT. ALL THE 
POSITIONS ARE REFERENCED TO THE ACOUSTIC CZNTER, WITH 
THE Z COMPONENT BEING MEASURED UPWARDS FROM THE ARRAY. 
DOUBLE PRECISION T(1000, i) (1000 3) gNT (10) 0) 
DOUBLE PRECISICN V,IC,D, £c NU MER, DENO 
D = 30.D0 
DO 11 I£ = 1,M 
NUMER = Ó. DO 
DENOM = 0. DO 
TC = T (1,4) 
56 02 J 1.3 
P [I J) = Ve v%(TC#TC-T (Ty J)**2) / (D*2. DO) 
NOMER = NUMER + P Í 
DENOM = DENOM + (ES EN "È 15. DO) #*2 
CONTINUE, 
NT (IL = Bd 
CONTINU I 
FETURN 
END 


105 


C 
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22 


SUEROUTINE POSN VC Pp 52) 


THIS FORTRAN SUBROUTINESI[I SEME NI SIE O IUS, I 
MZTHOD NAVY_A FOR ESE SMRPING APPARENT 2Csitiivire 


EYDROPRONES TIMES I TORSO OD 
IHE VELOCI TY V AT Tae 


APPARENT POSITION Pome se, 

Mm RAY @RANS ID Meio Nis ALLE 
ERENCED TO THE ACOUSTIC POENTON 
HE Z COMPONENT MEASURE UPWARDS 


INEUT ARE T 
LOCATIONS, 


DZVWOH rm | 
Metin ath 
KAH Ot3j t 
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L= = 

31 rrj G? 


GGA Otdgütidd 
0) 00 3 


OO EOOCHEG CH 
II 
ONNaWw HRM 
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iR 

bib HAOHO 
QOH DW try 
HH 

unu 

OO 

c Le 

«1i 

` = 

DI 

`a O 

ZO 

GO 

= 

by 

LETS 

- a 

Org 

t— 

= 3 

oo 
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. O 
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Qu 

CN Sel 

woe 

Q 
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OO UOH HAWAO 
HH 


O 
Ono 


=d (1) 
+ 


DO VEY HTC EE 
ME ( (I,J) A )) 


4 
C 
0 
O. 
2 
J P (I,J)*V* 

M = cg ; P (1,9 


J) 
ER = "véier" +P 
- TC*DSQRT(NUMEE/Di 


II AU O II ——Ə 


DS QRT (DCC) 


dch 


2 
OM) 


= *O 
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QO0O000001000 AQA 


ee Se eo SO ODI O OI 0.609 3919 19.49 £D O O JD O c C DO I 
SJERKROUITTNE POSLE Y; M, P, LSTI) 


DOOR ANS UDROUTINE IMPLEMENTS THE LEAST SQUARES 
POSTAR METAOD IS. 


INPUT ARE THE HYDROPHONE TIMES T FOR M SOUND SOURCE 
POSITIONS, AND THE VELOCITY V AT THE ARRAY. 
INPUTS AND OUTPUTS ARE THE SAME AS FOR THE 
SUBROUTINE POSNAV. 
INTEGER M,I,J 
DOUBLE PRÉCÍSION T(1000,4),P(1000,3),LST(1) 00) 
DOUBLE PRECISICN V,DISC, TC 
DO 11 I = 1,M 
ESE = 0-0 
Y S123) Zë 4) -T(1,3) 
DIS E = DTSC + br dy ee) 
22 CONTINUE 
EST (7) 9 (T (1,1) * (1,2) *T (1,3) - T (1,9) ) /2. DO 
Ee e 
P(I,J) = (V*LST(I)*P(I,J))/DSQRT (DISC 
3 coni nns! ( ADO ( ) 
1 CONTINUE 
FETURN 
END 
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SUBROUTINE POSISC (T,V,H,P,LST) 


THIS FORTRAN SUBROUTINE IMPLEMENTS THE BIAS ADJUSTED 
LEAST SOUARES ELANAR METHOD L.S.C. FOR ESTIMATING 
APPARENT POSITIONS. 
INPUTS AND OUTPUTS ARE THE SAME AS FOR THE 
SUBROUTINE POSIS. 
INTEGER I,M,J,K,L 
DOUBLE PRÉCÍSÍCÁÓ T(1000 UY 00, 3) ¿15 z (190) 
DOUBLE PRECISICN V M us DISC, TD,RÍSQR,R2SQR,E(3) 
SET UP HYDROPHCNE POSITIONS 
DO 44 J =1,3 
pO 55 I= 1,4 
A(T J) =°-15.D0 
CONTINUE 
A (J J) = 15.D0 
CONTINUE 
CAICULATE ORIGINAL L.S. SOLUTION 
DO 11 I= 1,5 
DISC 5807979 
e $3) lesa m ET 
DISC = DISC + Lu dt 
CONTINUE 
ck 3) A D DO 
zd # 
P(I,J) = (V*LST(I)*P(I,J)) /DSORT(DISC 
FORCE, — e 2 a Gg | 
= + + 
D (8) = 0:00 ' Í 
Se ATH z D à) (Psa A (4,K))**2 
= + = 
Maa a 
DESC = 0 D0 
TD = 0.D0 
DO 77 J = 1,3 
D (J) - 0:D0 
n SÉ š P3 (P(I,K A (J,K)) **2 
= + ue 
EE PH ies ` 
DISC = DISC 4 Wii -p (a )**2 
TD = TD + D(J) 
CONTINUE 
CALCULATE BIAS VECTOR E 
DISC = DSORT (DISC) 
SE ttp D (ü) ) *(D(8) -D(J) )/Z(DISCX*2. DO) ) -P (I,J) 
CONTINUE ELLA, í 
R2 SOR = 0.D0 
DO 99 J =1,3 
P(I,J) = SU me E 
a2SQR - R2SQR + P(1,3)**2 
CONTINUE 
ADJUST ORIGINAL RAY TRANSIT TIME 
LST (I) = LST (1) * (DS ORT (R2SQR/R1SQR)) 
CONTINUE 
RETURN 
END 
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SUBROUTINE POSMLP (T,V,M,P,MT,VAR) 


THIS FORTRAN SUBROUTINE IMPLEMENTS THE MAXC MUI 
LIKELIHOOD PLANAR METHOD M.-L. P. FOR ESTIMATING 
APPARENT POSITIONS. 
INPUTS ARE THE HYDROPHONE TIMES T FOR M SOURCE 
POSITIONS, AND THE VELOCITY OF SOUND V AT THE ARRAY. 
OUTPUTS ARE THE APPARENT POSITION ESTIMATES P AND 
RAY TRANSIT TIMES MT FOR EACH OF THE M POSITIONS. 
ALSO OUTPUT IS A VECTOR VAR OF M ESTIMATES OF THE 
TIMING ERROR VARIANCE. POSITIONS ARE ALL 2 EFERENCED 
TO THE ACOUSTIC CENTER OF THE ARRAY, WITH THE Z 
COMPONENT MEASURED UPWARDS FROM THE ARRAY. 
INTEGER M, I, J 
DOUBLE PRECISICN T (19. 00,4) ¢P( 1900, 3) 37 | 0) 
DOUBLE PRECISION V,TC,D,DIFF, NUMER, DENOM,T) L,C (3) 
DOUBLE PRECISICN VAR( 1000) ,cO (3) ,X(3) 
D = 30.D0 
INITIALIZE THE DIRECTION COSINE ESTIMATE Bi 
USING THE LEAST SQUARES SOLUTION. 
LOCP THROUGH THE M POSITIONS 
DO 11 T = 1,M 
= T(1,4) 
DENOM = (TC-1(1,1))**2+ (TC-T (L, 2) y ** 2 
TC-T(1,3)) **2 
ps E š It T(I,J))/DSQRT (DENOM) 
COSE us i ü 
ITER = 0 
DENOM = 0.DO. 
WERE ITER + 1 : 
TC = 0.D0 
USE ITERATICN FORMULAE TO DEVELOPE COSINES 
AND TIME VALUES. 
DO 33. J.= 1,3 
CO (J) = € (2) 
DENOM = DENOM * T(I,J)*C(J) 
TC = TC + T(I,J) 
CONTINUE 
TC = ((TC+T (1,4))+D*( ¿nic 42t 3)) 71) 74.00 
DENOM = DENCH - TO * (C (1) #C ( 
DO 66 j = 1 
Gig) = SEE, d: TC) /D ENOM 
mM = (J) C [g)) 
CONTINUE 
CHECK ITERATION TOLERANCES 
DIFF = DMAX1(X(1 x (2), x (3 
IF { ITER - it 106 Y), 
ORMAT Et d$ AP TONS EX CEED 719} 
FORMAT (2X,* IN POSITION NO. gé 
WRITE(6 60 IT 
WRITEIS.61 
G9 TO 11 
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aaa 


70 


DU 


55 
11 


IF ( TOL .Li. BDIPP ) Gomme 22 
NUMER = 0.D0 
DENOM = 0.DO 
CONVERT COSINES TO POSITIONS, AND TRANS. ATE TO 
ACCUSTIC CENTER REFERENCED COORDINATES. 
DO 44 J = 1,3 
P (I J) = ° W*C (J) *TC 
DENÓM = DENOM + P(I,J)**2 
P(I,J) = P(I,J) ~%5.D0 
NUMER = NUMER + P(I,J)**2 
CONTINUE 


MT (I) = TC*DSQRT (NUMER/DENOM) 
ESTIMATE VARTANCE 
VAR(I) - 0.DQ 
DO 55 J = 1,23 
VAR (I) = VAR(I) + (T (1,3) -TC+D*C(JI) /T)**2 
CONTINUE 
= VAR(I)/4. DO 
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AAD 


000 


DSDT = DSDIFI (I) # (DC (I) - V*TAU) /VK (I) 
L = di ( 11/0508 (K (I) ) 
S = S+TK(I)/DSORT(K(I)) 
11 CONTINUE 
DO =" s 
b d i J)s (TK (I) *DLDC (J) ) / LL*L*DSQRT(K (2) )) 
E 7 
33 cont INGE 
du Eszen (V*TC-D*L 
GP(4,1)=(GP (4,1) -D*DLDC (1) *(1.D0- SU. er 
GENT, Bio i a -DLDC zx ry AV EEK (3 
Ge (I, T 1. DO {GP AU (I) ¥L¥L)} 
GP = (T(L) *L* ( STA LDS 1)) 


I,4) = GP (1,4) Y (VK(I *L*L L es 
22 CONTINUE i c | 


GP (4,4 Jd DSDT4(W*TC-DHMeE))-—D*DDDOGSHMU ES 
(4,4) ( ( )) -D#DLDT# (1. 20-5) ) 


INVERT DERIVATIVE MATRIX 
CALL GAUSS3 (4,0,GP, GPi, IER, 4) 
CALCULATE INITIAL NEWTON SEARCH SOLUTION 
40 DO 44 I 
SCI ya 
44 CONTINU 


aaa 090 


etn- -TK (I) / (DSQRT(K(I)) *1) 
TAU-(V*TC-D*L)/(V* (1.D0-S)) 


U 

E GPI(4.I) *GC (I 
== = e He 
2c 


uu 
(nor 
OQ 
O 


4) - GPI(4,4) *GC (4) 
PREPARE FOR GOLDEN SECTION SEA RTII 


EPS = DM AX ID) SES AAA 
ia Le ) « (EPS/ )) 
R*(B(i)-A(I 
) E cd } + (B(1)-4 (1)) 
66 CONTINUE 
A 3) - TAU 
X (5) ) = O r) 


START o pete LOOP ( COMMENCE GOLDEN ae 


ON SEARCH 
PROVE THE INITIAL NEWTON SOLU N 


m 
ION ) 
TOO eO 
70 CALL OBJFCN JF14 I, K TK,5, X1 ,D, V, T, TC) 
+ 


50 ) GO TO 400 


410 FORM EE G.S Se T 5 DS rO) 


ht + 


G 
400 DIFE = 


) «8 tr} = ay 


È 
717 CONT 


X UM à (3) + B(4) - X71 (4) 
TAU =_X2 (4 
3 DIFF - DSQRT(DIFF + (A (4)-B(4))**2) 
c CHECK TOLERANCES - STOP G.S.S. IF TIGHT FNOJGH 
E IEPS ET. DIE) GO TO 100 
' CALL OBJFCN (F2,L,K, TK,S,X2,D,V,T,TC) 
C CHOOSE IMPROVING DIRECTION 
IE (F1 «LT. F2) 30 TO 90 
Hm. xdi 
C = X1(1) 
= Ia) = m 
x1 (4) = x 
b 297 0! 
C 
90 DO 99 I = 1,3 
spl = X (I 
XT[I) -,A (H +B(I)-X1(1) 
C(t) = xi 
99 CONTINUE 
41d) = Ad) +30) -x1(8) 
GO TO 70 
€ 
c CHECK TOLERANCES FOR NEWTON SEARCH ITERATID NS, 
C AND PREP FOR NEXT NEWTON ITERATION IF NECESSARY 
I. "H y TAs (c(t) C1(I)) 
LEE cu PPM 
IF ( TOL] Lf. pitt $ 0 10 
DIFF = DABS (TAU-T1 
È IF ( TOL2 Stats Go! TO 10 SI 
c DONE WITH  II-TH SET OF TIMES AND POSITIONS. MAKE 
C ESTIMATES, AND GC ON TO NEXT POSITION TO BE ESTIMATED 
VARSL(II) = o. pp 
O Ban 
NUMER = 0. DO 
DO 155 J = 1,3 
VARML (II) = VARUL(II) + TK(J)*TK (J) 
È ) = V*TAU*FC(J 
D yo - DENOM * P(II,J)**2 
P J) = VSTAU*C [3] 7 D72: DO 
NOLER = NUMER 4 Pl 6D 
155 CONTINUE 
VARHL(II) -. (VARML (II) fr: TAU) **2) /4.D0 
SDEV (LI T) = DSCRT (VA RML (21) 
ULT (II) = TALFDSORT (NUMER/DENON) 
5 CONTI 
C WRITE 7 4900) (SDEV(I),I = 1,5) 
C900 FORMA 6 E15. 5) 
END 
C 
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C 
C 
C 
C 
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€ 
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SUBROUTINE OBJ PCN (27 tye 75 AUN O) 


THIS FORTRAN SUBROUTINE IS CALLED BY THE SJ BROUTINE 
POSMLS. IT CALCULATES NEW VALUES FOR SEVERAL 
VARIABLES AS WELL AS A FUNCTIONAL VALUE WHICH IS THE 
DECISION FACTOR DETER MINING THE APPROPRIATE 
IMERCVING DIRECTION FOR THE GOLD2N SECTION SEARCH. 
INTEGER I 
DOUBLE PRECISION F,L,K (BLAKE (3) SUOLO) 
S = 0.D0 
L = 0.D0 
ni K (I) = {hy **2 (D/V) **2-((2. DO*D*X (4) *X(1))/V) 
= * - a 
TE LA) = uer 
S = S + a 1n E ( 5 ) 
L = L + (X(1)*Tx(1))/DSOAI (K (1)) 
CONTINUE 
F = 0.D0 
d £4 Te? (x TK (I) / (DSQRT (K (I) ) *L) ) ** 2 
= + == Š ` Ç L 

F S E+ (X (8) (V*TC-D*L)/(V#*(1.D0-S))) #*2 

= + [> = - ë - 
RETURN 
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